Branch data Line data Source code
1 : : /*
2 : : * spline.cpp - spline class implementation
3 : : *
4 : : * Copyright (C) 2005, 2006, 2009 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : : #include <assert.h>
33 : :
34 : : #include "logging.h"
35 : : #include "complex.h"
36 : : #include "object.h"
37 : : #include "vector.h"
38 : : #include "poly.h"
39 : : #include "tvector.h"
40 : : #include "tridiag.h"
41 : : #include "spline.h"
42 : :
43 : : namespace qucs {
44 : :
45 : : // Constructor creates an instance of the spline class.
46 : 0 : spline::spline () {
47 : 0 : x = f0 = f1 = f2 = f3 = NULL;
48 : 0 : d0 = dn = 0;
49 : 0 : n = 0;
50 : 0 : boundary = SPLINE_BC_NATURAL;
51 : 0 : }
52 : :
53 : : // Constructor creates an instance of the spline class with given boundary.
54 : 0 : spline::spline (int b) {
55 : 0 : x = f0 = f1 = f2 = f3 = NULL;
56 : 0 : d0 = dn = 0;
57 : 0 : n = 0;
58 : 0 : boundary = b;
59 : 0 : }
60 : :
61 : : // Constructor creates an instance of the spline class with vector data given.
62 : 0 : spline::spline (vector y, vector t) {
63 : 0 : x = f0 = f1 = f2 = f3 = NULL;
64 : 0 : d0 = dn = 0;
65 : 0 : n = 0;
66 : 0 : boundary = SPLINE_BC_NATURAL;
67 [ # # ][ # # ]: 0 : vectors (y, t);
[ # # ][ # # ]
[ # # ][ # # ]
68 : 0 : construct ();
69 : 0 : }
70 : :
71 : : // Constructor creates an instance of the spline class with tvector data given.
72 : 0 : spline::spline (tvector<nr_double_t> y, tvector<nr_double_t> t) {
73 : 0 : x = f0 = f1 = f2 = f3 = NULL;
74 : 0 : d0 = dn = 0;
75 : 0 : n = 0;
76 : 0 : boundary = SPLINE_BC_NATURAL;
77 [ # # ][ # # ]: 0 : vectors (y, t);
[ # # ][ # # ]
78 : 0 : construct ();
79 : 0 : }
80 : :
81 : : #define t_ (t)
82 : : #define y_ (y)
83 : :
84 : : // Pass interpolation datapoints as vectors.
85 : 0 : void spline::vectors (vector y, vector t) {
86 : 0 : int i = t.getSize ();
87 [ # # ][ # # ]: 0 : assert (y.getSize () == i && i >= 3);
[ # # ]
88 : :
89 : : // create local copy of f(x)
90 : 0 : realloc (i);
91 [ # # ]: 0 : for (i = 0; i <= n; i++) {
92 : 0 : f0[i] = real (y_(i)); x[i] = real (t_(i));
93 : : }
94 : 0 : }
95 : :
96 : : // Pass interpolation datapoints as tvectors.
97 : 0 : void spline::vectors (tvector<nr_double_t> y, tvector<nr_double_t> t) {
98 : 0 : int i = t.getSize ();
99 [ # # ][ # # ]: 0 : assert (y.getSize () == i && i >= 3);
[ # # ]
100 : :
101 : : // create local copy of f(x)
102 : 0 : realloc (i);
103 [ # # ]: 0 : for (i = 0; i <= n; i++) {
104 : 0 : f0[i] = y_(i); x[i] = t_(i);
105 : : }
106 : 0 : }
107 : :
108 : : // Pass interpolation datapoints as pointers.
109 : 0 : void spline::vectors (nr_double_t * y, nr_double_t * t, int len) {
110 : 0 : int i = len;
111 [ # # ]: 0 : assert (i >= 3);
112 : :
113 : : // create local copy of f(x)
114 : 0 : realloc (i);
115 [ # # ]: 0 : for (i = 0; i <= n; i++) {
116 : 0 : f0[i] = y[i]; x[i] = t[i];
117 : : }
118 : 0 : }
119 : :
120 : : // Reallocate vector data if necessary.
121 : 0 : void spline::realloc (int size) {
122 [ # # ]: 0 : if (n != size - 1) {
123 : 0 : n = size - 1;
124 [ # # ][ # # ]: 0 : if (f0) delete[] f0;
125 : 0 : f0 = new nr_double_t[n+1];
126 [ # # ][ # # ]: 0 : if (x) delete[] x;
127 : 0 : x = new nr_double_t[n+1];
128 : : }
129 [ # # ][ # # ]: 0 : if (f1) { delete[] f1; f1 = NULL; }
130 [ # # ][ # # ]: 0 : if (f2) { delete[] f2; f2 = NULL; }
131 [ # # ][ # # ]: 0 : if (f3) { delete[] f3; f3 = NULL; }
132 : 0 : }
133 : :
134 : : // Construct cubic spline interpolation coefficients.
135 : 0 : void spline::construct (void) {
136 : :
137 : : // calculate first derivative h = f'(x)
138 : : int i;
139 : 0 : nr_double_t * h = new nr_double_t[n+1];
140 [ # # ]: 0 : for (i = 0; i < n; i++) {
141 : 0 : h[i] = x[i+1] - x[i];
142 [ # # ]: 0 : if (h[i] == 0.0) {
143 : : logprint (LOG_ERROR, "ERROR: Duplicate points in spline: %g, %g\n",
144 : 0 : x[i], x[i+1]);
145 : : }
146 : : }
147 : :
148 : : // first kind of cubic splines
149 [ # # ][ # # ]: 0 : if (boundary == SPLINE_BC_NATURAL || boundary == SPLINE_BC_CLAMPED) {
150 : :
151 : : // setup right hand side
152 : 0 : nr_double_t * b = new nr_double_t[n+1]; // b as in Ax = b
153 [ # # ]: 0 : for (i = 1; i < n; i++) {
154 : 0 : nr_double_t _n = f0[i+1] * h[i-1] - f0[i] * (h[i] + h[i-1]) +
155 : 0 : f0[i-1] * h[i];
156 : 0 : nr_double_t _d = h[i-1] * h[i];
157 : 0 : b[i] = 3 * _n / _d;
158 : : }
159 [ # # ]: 0 : if (boundary == SPLINE_BC_NATURAL) {
160 : : // natural boundary condition
161 : 0 : b[0] = 0;
162 : 0 : b[n] = 0;
163 [ # # ]: 0 : } else if (boundary == SPLINE_BC_CLAMPED) {
164 : : // start and end derivatives given
165 : 0 : b[0] = 3 * ((f0[1] - f0[0]) / h[0] - d0);
166 : 0 : b[n] = 3 * (dn - (f0[n] - f0[n-1]) / h[n-1]);
167 : : }
168 : :
169 : 0 : nr_double_t * u = new nr_double_t[n+1];
170 : 0 : nr_double_t * z = b; // reuse storage
171 [ # # ]: 0 : if (boundary == SPLINE_BC_NATURAL) {
172 : 0 : u[0] = 0;
173 : 0 : z[0] = 0;
174 [ # # ]: 0 : } else if (boundary == SPLINE_BC_CLAMPED) {
175 : 0 : u[0] = h[0] / (2 * h[0]);
176 : 0 : z[0] = b[0] / (2 * h[0]);
177 : : }
178 : :
179 [ # # ]: 0 : for (i = 1; i < n; i++) {
180 : 0 : nr_double_t p = 2 * (h[i] + h[i-1]) - h[i-1] * u[i-1]; // pivot
181 : 0 : u[i] = h[i] / p;
182 : 0 : z[i] = (b[i] - z[i-1] * h[i-1]) / p;
183 : : }
184 [ # # ]: 0 : if (boundary == SPLINE_BC_NATURAL) {
185 : 0 : z[n] = 0;
186 [ # # ]: 0 : } else if (boundary == SPLINE_BC_CLAMPED) {
187 : 0 : nr_double_t p = h[n-1] * (2 - u[n-1]);
188 : 0 : z[n] = (b[n] - z[n-1] * h[n-1]) / p;
189 : : }
190 : :
191 : : // back substitution
192 : 0 : f1 = u; // reuse storage
193 : 0 : f2 = z;
194 : 0 : f3 = h;
195 : 0 : f2[n] = z[n];
196 : 0 : f3[n] = 0;
197 [ # # ]: 0 : for (i = n - 1; i >= 0; i--) {
198 : 0 : f2[i] = z[i] - u[i] * f2[i+1];
199 : 0 : f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (f2[i+1] + 2 * f2[i]) / 3;
200 : 0 : f3[i] = (f2[i+1] - f2[i]) / (3 * h[i]);
201 : : }
202 : :
203 : : // set up last slot for extrapolation above
204 [ # # ]: 0 : if (boundary == SPLINE_BC_NATURAL) {
205 : 0 : f1[n] = f1[n-1] + (x[n] - x[n-1]) * f2[n-1];
206 [ # # ]: 0 : } else if (boundary == SPLINE_BC_CLAMPED) {
207 : 0 : f1[n] = dn;
208 : : }
209 : 0 : f2[n] = 0;
210 : 0 : f3[n] = 0;
211 : : }
212 : :
213 : : // second kind of cubic splines
214 [ # # ]: 0 : else if (boundary == SPLINE_BC_PERIODIC) {
215 : : // non-trigdiagonal equations - periodic boundary condition
216 : 0 : nr_double_t * z = new nr_double_t[n+1];
217 [ # # ]: 0 : if (n == 2) {
218 : 0 : nr_double_t B = h[0] + h[1];
219 : 0 : nr_double_t A = 2 * B;
220 : : nr_double_t b[2], det;
221 : 0 : b[0] = 3 * ((f0[2] - f0[1]) / h[1] - (f0[1] - f0[0]) / h[0]);
222 : 0 : b[1] = 3 * ((f0[1] - f0[2]) / h[0] - (f0[2] - f0[1]) / h[1]);
223 : 0 : det = 3 * B * B;
224 : 0 : z[1] = ( A * b[0] - B * b[1]) / det;
225 : 0 : z[2] = (-B * b[0] + A * b[1]) / det;
226 : 0 : z[0] = z[2];
227 : : }
228 : : else {
229 : 0 : tridiag<nr_double_t> sys;
230 [ # # ]: 0 : tvector<nr_double_t> o (n);
231 [ # # ]: 0 : tvector<nr_double_t> d (n);
232 : 0 : tvector<nr_double_t> b;
233 : 0 : b.setData (&z[1], n);
234 [ # # ]: 0 : for (i = 0; i < n - 1; i++) {
235 [ # # ]: 0 : o(i) = h[i+1];
236 [ # # ]: 0 : d(i) = 2 * (h[i+1] + h[i]);
237 [ # # ]: 0 : b(i) = 3 * ((f0[i+2] - f0[i+1]) / h[i+1] - (f0[i+1] - f0[i]) / h[i]);
238 : : }
239 [ # # ]: 0 : o(i) = h[0];
240 [ # # ]: 0 : d(i) = 2 * (h[0] + h[i]);
241 [ # # ]: 0 : b(i) = 3 * ((f0[1] - f0[i+1]) / h[0] - (f0[i+1] - f0[i]) / h[i]);
242 : 0 : sys.setDiagonal (&d);
243 : 0 : sys.setOffDiagonal (&o);
244 : 0 : sys.setRHS (&b);
245 : 0 : sys.setType (TRIDIAG_SYM_CYCLIC);
246 [ # # ]: 0 : sys.solve ();
247 : 0 : z[0] = z[n];
248 : : }
249 : :
250 : 0 : f1 = new nr_double_t[n+1];
251 : 0 : f2 = z; // reuse storage
252 : 0 : f3 = h;
253 [ # # ]: 0 : for (i = n - 1; i >= 0; i--) {
254 : 0 : f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (z[i+1] + 2 * z[i]) / 3;
255 : 0 : f3[i] = (z[i+1] - z[i]) / (3 * h[i]);
256 : : }
257 : 0 : f1[n] = f1[0];
258 : 0 : f2[n] = f2[0];
259 : 0 : f3[n] = f3[0];
260 : : }
261 : 0 : }
262 : :
263 : : // Returns pointer to the first value greater than the given one.
264 : 0 : nr_double_t * spline::upper_bound (nr_double_t * first, nr_double_t * last,
265 : : nr_double_t value) {
266 : 0 : int half, len = last - first;
267 : : nr_double_t * centre;
268 : :
269 [ # # ]: 0 : while (len > 0) {
270 : 0 : half = len >> 1;
271 : 0 : centre = first;
272 : 0 : centre += half;
273 [ # # ]: 0 : if (value < *centre)
274 : 0 : len = half;
275 : : else {
276 : 0 : first = centre;
277 : 0 : ++first;
278 : 0 : len = len - half - 1;
279 : : }
280 : : }
281 : 0 : return first;
282 : : }
283 : :
284 : : // Evaluates the spline at the given position.
285 : 0 : poly spline::evaluate (nr_double_t t) {
286 : :
287 : : #ifndef PERIOD_DISABLED
288 [ # # ]: 0 : if (boundary == SPLINE_BC_PERIODIC) {
289 : : // extrapolation easy: periodically
290 : 0 : nr_double_t period = x[n] - x[0];
291 [ # # ]: 0 : while (t > x[n]) t -= period;
292 [ # # ]: 0 : while (t < x[0]) t += period;
293 : : }
294 : : #endif /* PERIOD_DISABLED */
295 : :
296 : 0 : nr_double_t * here = upper_bound (x, x+n+1, t);
297 : : nr_double_t y0, y1, y2;
298 [ # # ]: 0 : if (here == x) {
299 : 0 : nr_double_t dx = t - x[0];
300 : 0 : y0 = f0[0] + dx * f1[0];
301 : 0 : y1 = f1[0];
302 : 0 : return poly (t, y0, y1);
303 : : }
304 : : else {
305 : 0 : int i = here - x - 1;
306 : 0 : nr_double_t dx = t - x[i];
307 : : // value
308 : 0 : y0 = f0[i] + dx * (f1[i] + dx * (f2[i] + dx * f3[i]));
309 : : // first derivative
310 : 0 : y1 = f1[i] + dx * (2 * f2[i] + 3 * dx * f3[i]);
311 : : // second derivative
312 : 0 : y2 = 2 * f2[i] + 6 * dx * f3[i];
313 : 0 : return poly (t, y0, y1, y2);
314 : : }
315 : : }
316 : :
317 : : // Destructor deletes an instance of the spline class.
318 : 0 : spline::~spline () {
319 [ # # ][ # # ]: 0 : if (x) delete[] x;
[ # # ][ # # ]
320 [ # # ][ # # ]: 0 : if (f0) delete[] f0;
[ # # ][ # # ]
321 [ # # ][ # # ]: 0 : if (f1) delete[] f1;
[ # # ][ # # ]
322 [ # # ][ # # ]: 0 : if (f2) delete[] f2;
[ # # ][ # # ]
323 [ # # ][ # # ]: 0 : if (f3) delete[] f3;
[ # # ][ # # ]
324 : 0 : }
325 : :
326 : : } // namespace qucs
|
