Qucs-core  0.0.19
spline.cpp
Go to the documentation of this file.
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