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 1618 by gezelter, Mon Sep 12 17:09:26 2011 UTC

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

Diff Legend

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