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 1855 by gezelter, Tue Apr 2 18:31:51 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;
# 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]));
238 <  
225 >  dt = t - data_[j].first;
226 >  return data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));  
227   }
228  
229  
# Line 248 | Line 236 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
236    //   pair containing value of spline at t and first derivative at t
237  
238    if (!generated) generate();
251 –  RealType dt;
239    
240 <  if ( t < data.front().first || t > data.back().first ) {    
241 <    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 <  }
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  
264 –  int j;
265 –
246    if (isUniform) {    
247      
248 <    j = max(0, min(n-1, int((t - data[0].first) * dx)));  
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 ) {
255 >      if ( t < data_[i].first ) {
256          j = i-1;
257          break;
258        }      
# Line 281 | Line 261 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
261    
262    //  Evaluate the cubic polynomial.
263    
264 <  dt = t - data[j].first;
264 >  dt = t - data_[j].first;
265  
266 <  RealType yval = data[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
267 <  RealType dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
266 >  yval = data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j]));
267 >  dydx = b[j] + dt*(2.0 * c[j] + 3.0 * dt * d[j]);
268    
269    return make_pair(yval, dydx);
270   }

Diff Legend

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