ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 1536 by gezelter, Wed Jan 5 14:49:05 2011 UTC

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

Diff Legend

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