ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/math/CubicSpline.cpp
(Generate patch)

Comparing branches/development/src/math/CubicSpline.cpp (file contents):
Revision 1475 by gezelter, Wed Jul 21 17:51:22 2010 UTC vs.
Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

# Line 36 | Line 36
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "math/CubicSpline.hpp"
44   #include "utils/simError.h"
45   #include <cmath>
46 + #include <cstdio>
47   #include <algorithm>
46 – #include <iostream>
48  
49   using namespace OpenMD;
50   using namespace std;
51  
52 < CubicSpline::CubicSpline() : generated(false), isUniform(true) {}
52 > CubicSpline::CubicSpline() : generated(false), isUniform(true) {
53 >  data_.clear();
54 > }
55  
56 < void CubicSpline::addPoint(RealType xp, RealType yp) {
57 <  data.push_back(make_pair(xp, yp));
56 > void CubicSpline::addPoint(const RealType xp, const RealType yp) {
57 >  data_.push_back(make_pair(xp, yp));
58   }
59  
60   void CubicSpline::addPoints(const vector<RealType>& xps,
# Line 66 | Line 69 | void CubicSpline::addPoints(const vector<RealType>& xp
69    }
70  
71    for (int i = 0; i < xps.size(); i++)
72 <    data.push_back(make_pair(xps[i], yps[i]));
72 >    data_.push_back(make_pair(xps[i], yps[i]));
73   }
74  
75   void CubicSpline::generate() {
76    // Calculate coefficients defining a smooth cubic interpolatory spline.
77    //
78    // class values constructed:
79 <  //   n   = number of data points.
79 >  //   n   = number of data_ points.
80    //   x   = vector of independent variable values
81    //   y   = vector of dependent variable values
82    //   b   = vector of S'(x[i]) values.
83    //   c   = vector of S"(x[i])/2 values.
84    //   d   = vector of S'''(x[i]+)/6 values (i < n).
85    // Local variables:  
86 <  
86 >
87    RealType fp1, fpn, h, p;
88    
89    // make sure the sizes match
90    
91 <  n = data.size();  
89 <  x.resize(n);
90 <  y.resize(n);
91 >  n = data_.size();  
92    b.resize(n);
93    c.resize(n);
94    d.resize(n);
# Line 97 | Line 98 | void CubicSpline::generate() {
98    bool sorted = true;
99    
100    for (int i = 1; i < n; i++) {
101 <    if ( (data[i].first - data[i-1].first ) <= 0.0 ) sorted = false;
101 >    if ( (data_[i].first - data_[i-1].first ) <= 0.0 ) sorted = false;
102    }
103    
104    // sort if necessary
105    
106 <  if (!sorted) sort(data.begin(), data.end());  
106 >  if (!sorted) sort(data_.begin(), data_.end());  
107    
107 –  // Copy spline data out to separate arrays:
108 –  
109 –  for (int i = 0; i < n; i++) {
110 –    x[i] = data[i].first;
111 –    y[i] = data[i].second;
112 –  }
113 –  
108    // Calculate coefficients for the tridiagonal system: store
109    // sub-diagonal in B, diagonal in D, difference quotient in C.  
110    
111 <  b[0] = data[1].first - data[0].first;
112 <  c[0] = (data[1].second - data[0].second) / b[0];
111 >  b[0] = data_[1].first - data_[0].first;
112 >  c[0] = (data_[1].second - data_[0].second) / b[0];
113    
114    if (n == 2) {
115  
116      // Assume the derivatives at both endpoints are zero. Another
117      // assumption could be made to have a linear interpolant between
118      // the two points.  In that case, the b coefficients below would be
119 <    // (data[1].second - data[0].second) / (data[1].first - data[0].first)
119 >    // (data_[1].second - data_[0].second) / (data_[1].first - data_[0].first)
120      // and the c and d coefficients would both be zero.
121      b[0] = 0.0;
122 <    c[0] = -3.0 * pow((data[1].second - data[0].second) /
123 <                      (data[1].first-data[0].first), 2);
124 <    d[0] = -2.0 * pow((data[1].second - data[0].second) /
125 <                      (data[1].first-data[0].first), 3);
122 >    c[0] = -3.0 * pow((data_[1].second - data_[0].second) /
123 >                      (data_[1].first-data_[0].first), 2);
124 >    d[0] = -2.0 * pow((data_[1].second - data_[0].second) /
125 >                      (data_[1].first-data_[0].first), 3);
126      b[1] = b[0];
127      c[1] = 0.0;
128      d[1] = 0.0;
129 <    dx = 1.0 / (data[1].first - data[0].first);
129 >    dx = 1.0 / (data_[1].first - data_[0].first);
130      isUniform = true;
131      generated = true;
132      return;
# Line 141 | Line 135 | void CubicSpline::generate() {
135    d[0] = 2.0 * b[0];
136    
137    for (int i = 1; i < n-1; i++) {
138 <    b[i] = data[i+1].first - data[i].first;
138 >    b[i] = data_[i+1].first - data_[i].first;
139      if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
140 <    c[i] = (data[i+1].second - data[i].second) / b[i];
140 >    c[i] = (data_[i+1].second - data_[i].second) / b[i];
141      d[i] = 2.0 * (b[i] + b[i-1]);
142    }
143    
144    d[n-1] = 2.0 * b[n-2];
145    
146    // Calculate estimates for the end slopes using polynomials
147 <  // that interpolate the data nearest the end.
147 >  // that interpolate the data_ nearest the end.
148    
149    fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]);
150    if (n > 3) fp1 = fp1 + b[0]*((b[0] + b[1]) * (c[2] - c[1]) /
151                                 (b[1] + b[2]) -
152 <                               c[1] + c[0]) / (data[3].first - data[0].first);
152 >                               c[1] + c[0]) / (data_[3].first - data_[0].first);
153    
154    fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
155  
156    if (n > 3)  fpn = fpn + b[n-2] *
157      (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
158 <     (c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data[n-1].first - data[n-4].first);
158 >     (c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data_[n-1].first - data_[n-4].first);
159    
160    
161    // Calculate the right hand side and store it in C.
# Line 187 | Line 181 | void CubicSpline::generate() {
181    // Calculate the coefficients defining the spline.
182    
183    for (int i = 0; i < n-1; i++) {
184 <    h = data[i+1].first - data[i].first;
184 >    h = data_[i+1].first - data_[i].first;
185      d[i] = (c[i+1] - c[i]) / (3.0 * h);
186 <    b[i] = (data[i+1].second - data[i].second)/h - h * (c[i] + h * d[i]);
186 >    b[i] = (data_[i+1].second - data_[i].second)/h - h * (c[i] + h * d[i]);
187    }
188    
189    b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
190    
191 <  if (isUniform) dx = 1.0 / (data[1].first - data[0].first);
191 >  if (isUniform) dx = 1.0 / (data_[1].first - data_[0].first);
192    
193    generated = true;
194    return;
# Line 211 | Line 205 | RealType CubicSpline::getValueAt(RealType t) {
205    if (!generated) generate();
206    RealType dt;
207    
208 <  if ( t < data[0].first || t > data[n-1].first ) {    
208 >  if ( t < data_[0].first || t > data_[n-1].first ) {    
209      sprintf( painCave.errMsg,
210               "CubicSpline::getValueAt was passed a value outside the range of the spline!\n");
211      painCave.severity = OPENMD_ERROR;
# Line 226 | Line 220 | RealType CubicSpline::getValueAt(RealType t) {
220  
221    if (isUniform) {    
222      
223 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
223 >    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
224  
225    } else {
226  
227      j = n-1;
228      
229      for (int i = 0; i < n; i++) {
230 <      if ( t < data[i].first ) {
230 >      if ( t < data_[i].first ) {
231          j = i-1;
232          break;
233        }      
# Line 242 | Line 236 | RealType CubicSpline::getValueAt(RealType t) {
236    
237    //  Evaluate the cubic polynomial.
238    
239 <  dt = t - data[j].first;
240 <  return data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
239 >  dt = t - data_[j].first;
240 >  return data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
241    
242   }
243  
# Line 259 | Line 253 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
253    if (!generated) generate();
254    RealType dt;
255    
256 <  if ( t < data.front().first || t > data.back().first ) {    
256 >  if ( t < data_.front().first || t > data_.back().first ) {    
257      sprintf( painCave.errMsg,
258               "CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n");
259      painCave.severity = OPENMD_ERROR;
# Line 274 | Line 268 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
268  
269    if (isUniform) {    
270      
271 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
271 >    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
272  
273    } else {
274  
275      j = n-1;
276      
277      for (int i = 0; i < n; i++) {
278 <      if ( t < data[i].first ) {
278 >      if ( t < data_[i].first ) {
279          j = i-1;
280          break;
281        }      
# Line 290 | Line 284 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
284    
285    //  Evaluate the cubic polynomial.
286    
287 <  dt = t - data[j].first;
287 >  dt = t - data_[j].first;
288  
289 <  RealType yval = data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
289 >  RealType yval = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
290    RealType dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
291    
292    return make_pair(yval, dydx);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines