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 1618 by gezelter, Mon Sep 12 17:09:26 2011 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>
48  
# Line 59 | Line 60 | void CubicSpline::addPoints(const vector<RealType>& xp
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");
65 <    painCave.severity = OPENMD_ERROR;
66 <    painCave.isFatal = 1;
67 <    simError();    
68 <  }
69 <
70 <  for (int i = 0; i < xps.size(); 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  
# Line 151 | Line 146 | void CubicSpline::generate() {
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    
159 –  
155    // Calculate the right hand side and store it in C.
156    
157    c[n-1] = 3.0 * (fpn - c[n-2]);
# Line 193 | Line 188 | void CubicSpline::generate() {
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 202 | Line 197 | RealType CubicSpline::getValueAt(RealType t) {
197    //   value of spline at t.
198    
199    if (!generated) generate();
205 –  RealType dt;
200    
201 <  if ( t < data_[0].first || t > data_[n-1].first ) {    
202 <    sprintf( painCave.errMsg,
209 <             "CubicSpline::getValueAt was passed a value outside the range of the spline!\n");
210 <    painCave.severity = OPENMD_ERROR;
211 <    painCave.isFatal = 1;
212 <    simError();    
213 <  }
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  
218 –  int j;
219 –
207    if (isUniform) {    
208      
209      j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
# Line 236 | Line 223 | RealType CubicSpline::getValueAt(RealType t) {
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]));
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 250 | 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();
253 –  RealType dt;
278    
279 <  if ( t < data_.front().first || t > data_.back().first ) {    
280 <    sprintf( painCave.errMsg,
257 <             "CubicSpline::getValueAndDerivativeAt was passed a value outside the range of the spline!\n");
258 <    painCave.severity = OPENMD_ERROR;
259 <    painCave.isFatal = 1;
260 <    simError();    
261 <  }
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  
266 –  int j;
267 –
285    if (isUniform) {    
286      
287      j = max(0, min(n-1, int((t - data_[0].first) * dx)));  
# Line 285 | Line 302 | pair<RealType, RealType> CubicSpline::getValueAndDeriv
302    
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