| 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). |
| 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 |
|
|
| 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 |
|
|
| 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]); |
| 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))); |
| 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])); |
| 240 |
< |
|
| 226 |
> |
return data_[j].second + dt*(b[j] + dt*(c[j] + dt*d[j])); |
| 227 |
|
} |
| 228 |
|
|
| 229 |
|
|
| 236 |
|
// pair containing value of spline at t and first derivative at t |
| 237 |
|
|
| 238 |
|
if (!generated) generate(); |
| 253 |
– |
RealType dt; |
| 239 |
|
|
| 240 |
< |
if ( t < data_.front().first || t > data_.back().first ) { |
| 241 |
< |
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 |
< |
} |
| 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 |
|
|
| 266 |
– |
int j; |
| 267 |
– |
|
| 246 |
|
if (isUniform) { |
| 247 |
|
|
| 248 |
|
j = max(0, min(n-1, int((t - data_[0].first) * dx))); |
| 263 |
|
|
| 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 |
|
} |
