Qucs-core
0.0.19
|
00001 /* 00002 * spline.cpp - spline class implementation 00003 * 00004 * Copyright (C) 2005, 2006, 2009 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include <stdio.h> 00030 #include <stdlib.h> 00031 #include <string.h> 00032 #include <assert.h> 00033 #include <vector> 00034 00035 #include "logging.h" 00036 #include "complex.h" 00037 #include "object.h" 00038 #include "vector.h" 00039 #include "poly.h" 00040 #include "tvector.h" 00041 #include "tridiag.h" 00042 #include "spline.h" 00043 00044 namespace qucs { 00045 00046 // Constructor creates an instance of the spline class. 00047 spline::spline () { 00048 x = f0 = f1 = f2 = f3 = NULL; 00049 d0 = dn = 0; 00050 n = 0; 00051 boundary = SPLINE_BC_NATURAL; 00052 } 00053 00054 // Constructor creates an instance of the spline class with given boundary. 00055 spline::spline (int b) { 00056 x = f0 = f1 = f2 = f3 = NULL; 00057 d0 = dn = 0; 00058 n = 0; 00059 boundary = b; 00060 } 00061 00062 // Constructor creates an instance of the spline class with vector data given. 00063 spline::spline (qucs::vector y, qucs::vector t) { 00064 x = f0 = f1 = f2 = f3 = NULL; 00065 d0 = dn = 0; 00066 n = 0; 00067 boundary = SPLINE_BC_NATURAL; 00068 vectors (y, t); 00069 construct (); 00070 } 00071 00072 // Constructor creates an instance of the spline class with tvector data given. 00073 spline::spline (::std::vector<nr_double_t> y, ::std::vector<nr_double_t> t) { 00074 x = f0 = f1 = f2 = f3 = NULL; 00075 d0 = dn = 0; 00076 n = 0; 00077 boundary = SPLINE_BC_NATURAL; 00078 vectors (y, t); 00079 construct (); 00080 } 00081 00082 // Constructor creates an instance of the spline class with tvector data given. 00083 spline::spline (tvector<nr_double_t> y, tvector<nr_double_t> t) { 00084 x = f0 = f1 = f2 = f3 = NULL; 00085 d0 = dn = 0; 00086 n = 0; 00087 boundary = SPLINE_BC_NATURAL; 00088 vectors (y, t); 00089 construct (); 00090 } 00091 00092 #define t_ (t) 00093 #define y_ (y) 00094 00095 // Pass interpolation datapoints as vectors. 00096 void spline::vectors (qucs::vector y, qucs::vector t) { 00097 int i = t.getSize (); 00098 assert (y.getSize () == i && i >= 3); 00099 00100 // create local copy of f(x) 00101 realloc (i); 00102 for (i = 0; i <= n; i++) { 00103 f0[i] = real (y_(i)); x[i] = real (t_(i)); 00104 } 00105 } 00106 00107 // Pass interpolation datapoints as tvectors. 00108 void spline::vectors (::std::vector<nr_double_t> y, ::std::vector<nr_double_t> t) { 00109 int i = (int)t.size (); 00110 assert ((int)y.size () == i && i >= 3); 00111 00112 // create local copy of f(x) 00113 realloc (i); 00114 for (i = 0; i <= n; i++) { 00115 f0[i] = y[i]; x[i] = t[i]; 00116 } 00117 } 00118 00119 // Pass interpolation datapoints as tvectors. 00120 void spline::vectors (tvector<nr_double_t> y, tvector<nr_double_t> t) { 00121 int i = t.size (); 00122 assert (y.size () == i && i >= 3); 00123 00124 // create local copy of f(x) 00125 realloc (i); 00126 for (i = 0; i <= n; i++) { 00127 f0[i] = y_(i); x[i] = t_(i); 00128 } 00129 } 00130 00131 // Pass interpolation datapoints as pointers. 00132 void spline::vectors (nr_double_t * y, nr_double_t * t, int len) { 00133 int i = len; 00134 assert (i >= 3); 00135 00136 // create local copy of f(x) 00137 realloc (i); 00138 for (i = 0; i <= n; i++) { 00139 f0[i] = y[i]; x[i] = t[i]; 00140 } 00141 } 00142 00143 // Reallocate vector data if necessary. 00144 void spline::realloc (int size) { 00145 if (n != size - 1) { 00146 n = size - 1; 00147 delete[] f0; 00148 f0 = new nr_double_t[n+1]; 00149 delete[] x; 00150 x = new nr_double_t[n+1]; 00151 } 00152 delete[] f1; 00153 delete[] f2; 00154 delete[] f3; 00155 } 00156 00157 // Construct cubic spline interpolation coefficients. 00158 void spline::construct (void) { 00159 00160 // calculate first derivative h = f'(x) 00161 int i; 00162 nr_double_t * h = new nr_double_t[n+1]; 00163 for (i = 0; i < n; i++) { 00164 h[i] = x[i+1] - x[i]; 00165 if (h[i] == 0.0) { 00166 logprint (LOG_ERROR, "ERROR: Duplicate points in spline: %g, %g\n", 00167 x[i], x[i+1]); 00168 } 00169 } 00170 00171 // first kind of cubic splines 00172 if (boundary == SPLINE_BC_NATURAL || boundary == SPLINE_BC_CLAMPED) { 00173 00174 // setup right hand side 00175 nr_double_t * b = new nr_double_t[n+1]; // b as in Ax = b 00176 for (i = 1; i < n; i++) { 00177 nr_double_t _n = f0[i+1] * h[i-1] - f0[i] * (h[i] + h[i-1]) + 00178 f0[i-1] * h[i]; 00179 nr_double_t _d = h[i-1] * h[i]; 00180 b[i] = 3 * _n / _d; 00181 } 00182 if (boundary == SPLINE_BC_NATURAL) { 00183 // natural boundary condition 00184 b[0] = 0; 00185 b[n] = 0; 00186 } else if (boundary == SPLINE_BC_CLAMPED) { 00187 // start and end derivatives given 00188 b[0] = 3 * ((f0[1] - f0[0]) / h[0] - d0); 00189 b[n] = 3 * (dn - (f0[n] - f0[n-1]) / h[n-1]); 00190 } 00191 00192 nr_double_t * u = new nr_double_t[n+1]; 00193 nr_double_t * z = b; // reuse storage 00194 if (boundary == SPLINE_BC_NATURAL) { 00195 u[0] = 0; 00196 z[0] = 0; 00197 } else if (boundary == SPLINE_BC_CLAMPED) { 00198 u[0] = h[0] / (2 * h[0]); 00199 z[0] = b[0] / (2 * h[0]); 00200 } 00201 00202 for (i = 1; i < n; i++) { 00203 nr_double_t p = 2 * (h[i] + h[i-1]) - h[i-1] * u[i-1]; // pivot 00204 u[i] = h[i] / p; 00205 z[i] = (b[i] - z[i-1] * h[i-1]) / p; 00206 } 00207 if (boundary == SPLINE_BC_NATURAL) { 00208 z[n] = 0; 00209 } else if (boundary == SPLINE_BC_CLAMPED) { 00210 nr_double_t p = h[n-1] * (2 - u[n-1]); 00211 z[n] = (b[n] - z[n-1] * h[n-1]) / p; 00212 } 00213 00214 // back substitution 00215 f1 = u; // reuse storage 00216 f2 = z; 00217 f3 = h; 00218 f2[n] = z[n]; 00219 f3[n] = 0; 00220 for (i = n - 1; i >= 0; i--) { 00221 f2[i] = z[i] - u[i] * f2[i+1]; 00222 f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (f2[i+1] + 2 * f2[i]) / 3; 00223 f3[i] = (f2[i+1] - f2[i]) / (3 * h[i]); 00224 } 00225 00226 // set up last slot for extrapolation above 00227 if (boundary == SPLINE_BC_NATURAL) { 00228 f1[n] = f1[n-1] + (x[n] - x[n-1]) * f2[n-1]; 00229 } else if (boundary == SPLINE_BC_CLAMPED) { 00230 f1[n] = dn; 00231 } 00232 f2[n] = 0; 00233 f3[n] = 0; 00234 } 00235 00236 // second kind of cubic splines 00237 else if (boundary == SPLINE_BC_PERIODIC) { 00238 // non-trigdiagonal equations - periodic boundary condition 00239 ::std::vector<nr_double_t> z (n+1); 00240 if (n == 2) { 00241 nr_double_t B = h[0] + h[1]; 00242 nr_double_t A = 2 * B; 00243 nr_double_t b[2], det; 00244 b[0] = 3 * ((f0[2] - f0[1]) / h[1] - (f0[1] - f0[0]) / h[0]); 00245 b[1] = 3 * ((f0[1] - f0[2]) / h[0] - (f0[2] - f0[1]) / h[1]); 00246 det = 3 * B * B; 00247 z[1] = ( A * b[0] - B * b[1]) / det; 00248 z[2] = (-B * b[0] + A * b[1]) / det; 00249 z[0] = z[2]; 00250 } 00251 else { 00252 tridiag<nr_double_t> sys; 00253 ::std::vector<nr_double_t> o (n); 00254 ::std::vector<nr_double_t> d (n); 00255 ::std::vector<nr_double_t> b(&z[1],&z[n]); 00256 //b.setData (&z[1], n); 00257 for (i = 0; i < n - 1; i++) { 00258 o[i] = h[i+1]; 00259 d[i] = 2 * (h[i+1] + h[i]); 00260 b[i] = 3 * ((f0[i+2] - f0[i+1]) / h[i+1] - (f0[i+1] - f0[i]) / h[i]); 00261 z[i+1] = b[i]; 00262 } 00263 o[i] = h[0]; 00264 d[i] = 2 * (h[0] + h[i]); 00265 b[i] = 3 * ((f0[1] - f0[i+1]) / h[0] - (f0[i+1] - f0[i]) / h[i]); 00266 z[i+1] = b[i]; 00267 sys.setDiagonal (&d); 00268 sys.setOffDiagonal (&o); 00269 sys.setRHS (&b); 00270 sys.setType (TRIDIAG_SYM_CYCLIC); 00271 sys.solve (); 00272 z[0] = z[n]; 00273 } 00274 00275 f1 = new nr_double_t[n+1]; 00276 f2 = &z.front (); // reuse storage 00277 f3 = h; 00278 for (i = n - 1; i >= 0; i--) { 00279 f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (z[i+1] + 2 * z[i]) / 3; 00280 f3[i] = (z[i+1] - z[i]) / (3 * h[i]); 00281 } 00282 f1[n] = f1[0]; 00283 f2[n] = f2[0]; 00284 f3[n] = f3[0]; 00285 } 00286 } 00287 00288 // Returns pointer to the first value greater than the given one. 00289 nr_double_t * spline::upper_bound (nr_double_t * first, nr_double_t * last, 00290 nr_double_t value) { 00291 int half, len = last - first; 00292 nr_double_t * centre; 00293 00294 while (len > 0) { 00295 half = len >> 1; 00296 centre = first; 00297 centre += half; 00298 if (value < *centre) 00299 len = half; 00300 else { 00301 first = centre; 00302 ++first; 00303 len = len - half - 1; 00304 } 00305 } 00306 return first; 00307 } 00308 00309 // Evaluates the spline at the given position. 00310 poly spline::evaluate (nr_double_t t) { 00311 00312 #ifndef PERIOD_DISABLED 00313 if (boundary == SPLINE_BC_PERIODIC) { 00314 // extrapolation easy: periodically 00315 nr_double_t period = x[n] - x[0]; 00316 while (t > x[n]) t -= period; 00317 while (t < x[0]) t += period; 00318 } 00319 #endif /* PERIOD_DISABLED */ 00320 00321 nr_double_t * here = upper_bound (x, x+n+1, t); 00322 nr_double_t y0, y1, y2; 00323 if (here == x) { 00324 nr_double_t dx = t - x[0]; 00325 y0 = f0[0] + dx * f1[0]; 00326 y1 = f1[0]; 00327 return poly (t, y0, y1); 00328 } 00329 else { 00330 int i = here - x - 1; 00331 nr_double_t dx = t - x[i]; 00332 // value 00333 y0 = f0[i] + dx * (f1[i] + dx * (f2[i] + dx * f3[i])); 00334 // first derivative 00335 y1 = f1[i] + dx * (2 * f2[i] + 3 * dx * f3[i]); 00336 // second derivative 00337 y2 = 2 * f2[i] + 6 * dx * f3[i]; 00338 return poly (t, y0, y1, y2); 00339 } 00340 } 00341 00342 // Destructor deletes an instance of the spline class. 00343 spline::~spline () { 00344 if (x) delete[] x; 00345 delete[] f0; 00346 delete[] f1; 00347 delete[] f2; 00348 delete[] f3; 00349 } 00350 00351 } // namespace qucs