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 1479 by gezelter, Mon Jul 26 19:00:48 2010 UTC vs.
Revision 1877 by gezelter, Thu Jun 6 15:43:35 2013 UTC

# Line 35 | Line 35
35   *                                                                      
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).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
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"
43 – #include "utils/simError.h"
44   #include <cmath>
45 + #include <cassert>
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,
61                              const vector<RealType>& yps) {
62    
63 <  if (xps.size() != yps.size()) {
64 <    printf( painCave.errMsg,
65 <            "CubicSpline::addPoints was passed vectors of different length!\n");
66 <    painCave.severity = OPENMD_ERROR;
64 <    painCave.isFatal = 1;
65 <    simError();    
66 <  }
67 <
68 <  for (int i = 0; i < xps.size(); i++)
69 <    data.push_back(make_pair(xps[i], yps[i]));
63 >  assert(xps.size() == yps.size());
64 >  
65 >  for (unsigned int i = 0; i < xps.size(); i++)
66 >    data_.push_back(make_pair(xps[i], yps[i]));
67   }
68  
69   void CubicSpline::generate() {
70    // Calculate coefficients defining a smooth cubic interpolatory spline.
71    //
72    // class values constructed:
73 <  //   n   = number of data points.
73 >  //   n   = number of data_ points.
74    //   x   = vector of independent variable values
75    //   y   = vector of dependent variable values
76    //   b   = vector of S'(x[i]) values.
77    //   c   = vector of S"(x[i])/2 values.
78    //   d   = vector of S'''(x[i]+)/6 values (i < n).
79    // Local variables:  
80 <  
80 >
81    RealType fp1, fpn, h, p;
82    
83    // make sure the sizes match
84    
85 <  n = data.size();  
85 >  n = data_.size();  
86    b.resize(n);
87    c.resize(n);
88    d.resize(n);
# Line 95 | Line 92 | void CubicSpline::generate() {
92    bool sorted = true;
93    
94    for (int i = 1; i < n; i++) {
95 <    if ( (data[i].first - data[i-1].first ) <= 0.0 ) sorted = false;
95 >    if ( (data_[i].first - data_[i-1].first ) <= 0.0 ) sorted = false;
96    }
97    
98    // sort if necessary
99    
100 <  if (!sorted) sort(data.begin(), data.end());  
100 >  if (!sorted) sort(data_.begin(), data_.end());  
101    
102    // Calculate coefficients for the tridiagonal system: store
103    // sub-diagonal in B, diagonal in D, difference quotient in C.  
104    
105 <  b[0] = data[1].first - data[0].first;
106 <  c[0] = (data[1].second - data[0].second) / b[0];
105 >  b[0] = data_[1].first - data_[0].first;
106 >  c[0] = (data_[1].second - data_[0].second) / b[0];
107    
108    if (n == 2) {
109  
110      // Assume the derivatives at both endpoints are zero. Another
111      // assumption could be made to have a linear interpolant between
112      // the two points.  In that case, the b coefficients below would be
113 <    // (data[1].second - data[0].second) / (data[1].first - data[0].first)
113 >    // (data_[1].second - data_[0].second) / (data_[1].first - data_[0].first)
114      // and the c and d coefficients would both be zero.
115      b[0] = 0.0;
116 <    c[0] = -3.0 * pow((data[1].second - data[0].second) /
117 <                      (data[1].first-data[0].first), 2);
118 <    d[0] = -2.0 * pow((data[1].second - data[0].second) /
119 <                      (data[1].first-data[0].first), 3);
116 >    c[0] = -3.0 * pow((data_[1].second - data_[0].second) /
117 >                      (data_[1].first-data_[0].first), 2);
118 >    d[0] = -2.0 * pow((data_[1].second - data_[0].second) /
119 >                      (data_[1].first-data_[0].first), 3);
120      b[1] = b[0];
121      c[1] = 0.0;
122      d[1] = 0.0;
123 <    dx = 1.0 / (data[1].first - data[0].first);
123 >    dx = 1.0 / (data_[1].first - data_[0].first);
124      isUniform = true;
125      generated = true;
126      return;
# Line 132 | Line 129 | void CubicSpline::generate() {
129    d[0] = 2.0 * b[0];
130    
131    for (int i = 1; i < n-1; i++) {
132 <    b[i] = data[i+1].first - data[i].first;
132 >    b[i] = data_[i+1].first - data_[i].first;
133      if ( fabs( b[i] - b[0] ) / b[0] > 1.0e-5) isUniform = false;
134 <    c[i] = (data[i+1].second - data[i].second) / b[i];
134 >    c[i] = (data_[i+1].second - data_[i].second) / b[i];
135      d[i] = 2.0 * (b[i] + b[i-1]);
136    }
137    
138    d[n-1] = 2.0 * b[n-2];
139    
140    // Calculate estimates for the end slopes using polynomials
141 <  // that interpolate the data nearest the end.
141 >  // that interpolate the data_ nearest the end.
142    
143    fp1 = c[0] - b[0]*(c[1] - c[0])/(b[0] + b[1]);
144    if (n > 3) fp1 = fp1 + b[0]*((b[0] + b[1]) * (c[2] - c[1]) /
145                                 (b[1] + b[2]) -
146 <                               c[1] + c[0]) / (data[3].first - data[0].first);
146 >                               c[1] + c[0]) / (data_[3].first - data_[0].first);
147    
148    fpn = c[n-2] + b[n-2]*(c[n-2] - c[n-3])/(b[n-3] + b[n-2]);
149 <
149 >  
150    if (n > 3)  fpn = fpn + b[n-2] *
151 <    (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
152 <     (c[n-3] - c[n-4])/(b[n-3] + b[n-4]))/(data[n-1].first - data[n-4].first);
151 >                (c[n-2] - c[n-3] - (b[n-3] + b[n-2]) *
152 >                 (c[n-3] - c[n-4])/(b[n-3] + b[n-4])) /
153 >                (data_[n-1].first - data_[n-4].first);
154    
157 –  
155    // Calculate the right hand side and store it in C.
156    
157    c[n-1] = 3.0 * (fpn - c[n-2]);
# Line 178 | Line 175 | void CubicSpline::generate() {
175    // Calculate the coefficients defining the spline.
176    
177    for (int i = 0; i < n-1; i++) {
178 <    h = data[i+1].first - data[i].first;
178 >    h = data_[i+1].first - data_[i].first;
179      d[i] = (c[i+1] - c[i]) / (3.0 * h);
180 <    b[i] = (data[i+1].second - data[i].second)/h - h * (c[i] + h * d[i]);
180 >    b[i] = (data_[i+1].second - data_[i].second)/h - h * (c[i] + h * d[i]);
181    }
182    
183    b[n-1] = b[n-2] + h * (2.0 * c[n-2] + h * 3.0 * d[n-2]);
184    
185 <  if (isUniform) dx = 1.0 / (data[1].first - data[0].first);
185 >  if (isUniform) dx = 1.0 / (data_[1].first - data_[0].first);
186    
187    generated = true;
188    return;
189   }
190  
191 < RealType CubicSpline::getValueAt(RealType t) {
191 > RealType CubicSpline::getValueAt(const RealType& t) {
192    // Evaluate the spline at t using coefficients
193    //
194    // Input parameters
# Line 200 | Line 197 | RealType CubicSpline::getValueAt(RealType t) {
197    //   value of spline at t.
198    
199    if (!generated) generate();
203 –  RealType dt;
200    
201 <  if ( t < data[0].first || t > data[n-1].first ) {    
202 <    sprintf( painCave.errMsg,
207 <             "CubicSpline::getValueAt was passed a value outside the range of the spline!\n");
208 <    painCave.severity = OPENMD_ERROR;
209 <    painCave.isFatal = 1;
210 <    simError();    
211 <  }
201 >  assert(t > data_.front().first);
202 >  assert(t < data_.back().first);
203  
204    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
205    //  to t.
206  
216 –  int j;
217 –
207    if (isUniform) {    
208      
209 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
209 >    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
210  
211    } else {
212  
213      j = n-1;
214      
215      for (int i = 0; i < n; i++) {
216 <      if ( t < data[i].first ) {
216 >      if ( t < data_[i].first ) {
217          j = i-1;
218          break;
219        }      
# Line 233 | Line 222 | RealType CubicSpline::getValueAt(RealType t) {
222    
223    //  Evaluate the cubic polynomial.
224    
225 <  dt = t - data[j].first;
226 <  return data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
225 >  dt = t - data_[j].first;
226 >  return data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));  
227 > }
228 >
229 >
230 > void CubicSpline::getValueAt(const RealType& t, RealType& v) {
231 >  // Evaluate the spline at t using coefficients
232 >  //
233 >  // Input parameters
234 >  //   t = point where spline is to be evaluated.
235 >  // Output:
236 >  //   value of spline at t.
237    
238 +  if (!generated) generate();
239 +  
240 +  assert(t > data_.front().first);
241 +  assert(t < data_.back().first);
242 +
243 +  //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
244 +  //  to t.
245 +
246 +  if (isUniform) {    
247 +    
248 +    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
249 +
250 +  } else {
251 +
252 +    j = n-1;
253 +    
254 +    for (int i = 0; i < n; i++) {
255 +      if ( t < data_[i].first ) {
256 +        j = i-1;
257 +        break;
258 +      }      
259 +    }
260 +  }
261 +  
262 +  //  Evaluate the cubic polynomial.
263 +  
264 +  dt = t - data_[j].first;
265 +  v = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));  
266   }
267  
268  
269 < pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(RealType t) {
269 > pair<RealType, RealType> CubicSpline::getValueAndDerivativeAt(const RealType& t){
270    // Evaluate the spline and first derivative at t using coefficients
271    //
272    // Input parameters
# Line 248 | Line 275 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
275    //   pair containing value of spline at t and first derivative at t
276  
277    if (!generated) generate();
251 –  RealType dt;
278    
279 <  if ( t < data.front().first || t > data.back().first ) {    
280 <    sprintf( painCave.errMsg,
255 <             "CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n");
256 <    painCave.severity = OPENMD_ERROR;
257 <    painCave.isFatal = 1;
258 <    simError();    
259 <  }
279 >  assert(t > data_.front().first);
280 >  assert(t < data_.back().first);
281  
282    //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
283    //  to t.
284  
264 –  int j;
265 –
285    if (isUniform) {    
286      
287 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
287 >    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
288  
289    } else {
290  
291      j = n-1;
292      
293      for (int i = 0; i < n; i++) {
294 <      if ( t < data[i].first ) {
294 >      if ( t < data_[i].first ) {
295          j = i-1;
296          break;
297        }      
# Line 281 | Line 300 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
300    
301    //  Evaluate the cubic polynomial.
302    
303 <  dt = t - data[j].first;
303 >  dt = t - data_[j].first;
304  
305 <  RealType yval = data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
306 <  RealType dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
307 <  
305 >  yval = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
306 >  dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
307 >
308    return make_pair(yval, dydx);
309   }
310 +
311 + void CubicSpline::getValueAndDerivativeAt(const RealType& t, RealType& v,
312 +                                          RealType &dv) {
313 +  // Evaluate the spline and first derivative at t using coefficients
314 +  //
315 +  // Input parameters
316 +  //   t = point where spline is to be evaluated.
317 +
318 +  if (!generated) generate();
319 +  
320 +  assert(t > data_.front().first);
321 +  assert(t < data_.back().first);
322 +
323 +  //  Find the interval ( x[j], x[j+1] ) that contains or is nearest
324 +  //  to t.
325 +
326 +  if (isUniform) {    
327 +    
328 +    j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
329 +
330 +  } else {
331 +
332 +    j = n-1;
333 +    
334 +    for (int i = 0; i < n; i++) {
335 +      if ( t < data_[i].first ) {
336 +        j = i-1;
337 +        break;
338 +      }      
339 +    }
340 +  }
341 +  
342 +  //  Evaluate the cubic polynomial.
343 +  
344 +  dt = t - data_[j].first;
345 +
346 +  v = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
347 +  dv = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
348 + }

Diff Legend

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