# | 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); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |