Main Page | Class Hierarchy | Alphabetical List | Compound List | File List | Compound Members | File Members | Related Pages

qwt_spline.cpp

Go to the documentation of this file.
00001 /* -*- mode: C++ ; c-file-style: "stroustrup" -*- *****************************
00002  * Qwt Widget Library
00003  * Copyright (C) 1997   Josef Wilgen
00004  * Copyright (C) 2002   Uwe Rathmann
00005  * 
00006  * This library is free software; you can redistribute it and/or
00007  * modify it under the terms of the Qwt License, Version 1.0
00008  *****************************************************************************/
00009 
00010 #include "qwt_spline.h"
00011 #include "qwt_math.h"
00012 #include "qwt.h"
00013 
00018 double QwtSpline::value(double x) const
00019 {
00020     if (!d_a)
00021         return 0.0;
00022 
00023     const int i = lookup(x);
00024 
00025     const double delta = x - d_x[i];
00026     return( ( ( ( d_a[i] * delta) + d_b[i] ) 
00027         * delta + d_c[i] ) * delta + d_y[i] );
00028 }
00029 
00031 QwtSpline::QwtSpline()
00032 {
00033     d_a = d_b = d_c = 0;
00034     d_xbuffer = d_ybuffer = d_x = d_y = 0;
00035     d_size = 0;
00036     d_buffered = 0;
00037 }
00038 
00061 void QwtSpline::copyValues(int tf)
00062 {
00063     cleanup();
00064     d_buffered = tf;
00065 }
00066 
00068 QwtSpline::~QwtSpline()
00069 {
00070     cleanup();
00071 }
00072 
00074 int QwtSpline::lookup(double x) const
00075 {
00076     int i1, i2, i3;
00077     
00078     if (x <= d_x[0])
00079        i1 = 0;
00080     else if (x >= d_x[d_size - 2])
00081        i1 = d_size -2;
00082     else
00083     {
00084         i1 = 0;
00085         i2 = d_size -2;
00086         i3 = 0;
00087 
00088         while ( i2 - i1 > 1 )
00089         {
00090             i3 = i1 + ((i2 - i1) >> 1);
00091 
00092             if (d_x[i3] > x)
00093                i2 = i3;
00094             else
00095                i1 = i3;
00096 
00097         }
00098     }
00099     return i1;
00100 }
00101 
00102 
00127 int QwtSpline::recalc(double *x, double *y, int n, int periodic)
00128 {
00129     int i, rv = 0;
00130 
00131     cleanup();
00132 
00133     if (n > 2)
00134     {
00135         d_size = n;
00136 
00137         if (d_buffered)
00138         {
00139             d_xbuffer = new double[n-1];
00140             d_ybuffer = new double[n-1];
00141 
00142             if ((!d_xbuffer) || (!d_ybuffer))
00143             {
00144                 cleanup();
00145                 return Qwt::ErrNoMem;
00146             }
00147             else
00148             {
00149                 for (i=0;i<n;i++)
00150                 {
00151                     d_xbuffer[i] = x[i];
00152                     d_ybuffer[i] = y[i];
00153                 }
00154                 d_x = d_xbuffer;
00155                 d_y = d_ybuffer;
00156             }
00157         }
00158         else
00159         {
00160             d_x = x;
00161             d_y = y;
00162         }
00163         
00164         d_a = new double[n-1];
00165         d_b = new double[n-1];
00166         d_c = new double[n-1];
00167 
00168         if ( (!d_a) || (!d_b) || (!d_c) )
00169         {
00170             cleanup();
00171             return Qwt::ErrMono;
00172         }
00173 
00174         if(periodic)
00175            rv =  buildPerSpline();
00176         else
00177            rv =  buildNatSpline();
00178 
00179         if (rv) cleanup();
00180     }
00181 
00182     return rv;
00183 }
00184 
00194 int QwtSpline::buildNatSpline()
00195 {
00196     int i;
00197     double dy1, dy2;
00198     
00199     double *d = new double[d_size-1];
00200     double *h = new double[d_size-1];
00201     double *s = new double[d_size];
00202 
00203     if ( (!d) || (!h) || (!s) )
00204     {
00205         cleanup();
00206         if (h) delete[] h;
00207         if (s) delete[] s;
00208         if (d) delete[] d;
00209         return Qwt::ErrNoMem;
00210     }
00211 
00212     //
00213     //  set up tridiagonal equation system; use coefficient
00214     //  vectors as temporary buffers
00215     for (i=0; i<d_size - 1; i++) 
00216     {
00217         h[i] = d_x[i+1] - d_x[i];
00218         if (h[i] <= 0)
00219         {
00220             delete[] h;
00221             delete[] s;
00222             delete[] d;
00223             return Qwt::ErrMono;
00224         }
00225     }
00226     
00227     dy1 = (d_y[1] - d_y[0]) / h[0];
00228     for (i = 1; i < d_size - 1; i++)
00229     {
00230         d_b[i] = d_c[i] = h[i];
00231         d_a[i] = 2.0 * (h[i-1] + h[i]);
00232 
00233         dy2 = (d_y[i+1] - d_y[i]) / h[i];
00234         d[i] = 6.0 * ( dy1 - dy2);
00235         dy1 = dy2;
00236     }
00237 
00238     //
00239     // solve it
00240     //
00241     
00242     // L-U Factorization
00243     for(i = 1; i < d_size - 2;i++)
00244     {
00245         d_c[i] /= d_a[i];
00246         d_a[i+1] -= d_b[i] * d_c[i]; 
00247     }
00248 
00249     // forward elimination
00250     s[1] = d[1];
00251     for(i=2;i<d_size - 1;i++)
00252        s[i] = d[i] - d_c[i-1] * s[i-1];
00253     
00254     // backward elimination
00255     s[d_size - 2] = - s[d_size - 2] / d_a[d_size - 2];
00256     for (i= d_size -3; i > 0; i--)
00257        s[i] = - (s[i] + d_b[i] * s[i+1]) / d_a[i];
00258 
00259     //
00260     // Finally, determine the spline coefficients
00261     //
00262     s[d_size - 1] = s[0] = 0.0;
00263     for (i = 0; i < d_size - 1; i++)
00264     {
00265         d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00266         d_b[i] = 0.5 * s[i];
00267         d_c[i] = ( d_y[i+1] - d_y[i] ) 
00268             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00269     }
00270 
00271     delete[] d;
00272     delete[] s;
00273     delete[] h;
00274 
00275     return 0;
00276     
00277 }
00278 
00288 int QwtSpline::buildPerSpline()
00289 {
00290     int i,imax;
00291     double sum;
00292     double dy1, dy2,htmp;
00293     
00294     double *d = new double[d_size-1];
00295     double *h = new double[d_size-1];
00296     double *s = new double[d_size];
00297     
00298     if ( (!d) || (!h) || (!s) )
00299     {
00300         cleanup();
00301         if (h) delete[] h;
00302         if (s) delete[] s;
00303         if (d) delete[] d;
00304         return Qwt::ErrNoMem;
00305     }
00306 
00307     //
00308     //  setup equation system; use coefficient
00309     //  vectors as temporary buffers
00310     //
00311     for (i=0; i<d_size - 1; i++)
00312     {
00313         h[i] = d_x[i+1] - d_x[i];
00314         if (h[i] <= 0.0)
00315         {
00316             delete[] h;
00317             delete[] s;
00318             delete[] d;
00319             return Qwt::ErrMono;
00320         }
00321     }
00322     
00323     imax = d_size - 2;
00324     htmp = h[imax];
00325     dy1 = (d_y[0] - d_y[imax]) / htmp;
00326     for (i=0; i <= imax; i++)
00327     {
00328         d_b[i] = d_c[i] = h[i];
00329         d_a[i] = 2.0 * (htmp + h[i]);
00330         dy2 = (d_y[i+1] - d_y[i]) / h[i];
00331         d[i] = 6.0 * ( dy1 - dy2);
00332         dy1 = dy2;
00333         htmp = h[i];
00334     }
00335 
00336     //
00337     // solve it
00338     //
00339     
00340     // L-U Factorization
00341     d_a[0] = sqrt(d_a[0]);
00342     d_c[0] = h[imax] / d_a[0];
00343     sum = 0;
00344 
00345     for(i=0;i<imax-1;i++)
00346     {
00347         d_b[i] /= d_a[i];
00348         if (i > 0)
00349            d_c[i] = - d_c[i-1] * d_b[i-1] / d_a[i];
00350         d_a[i+1] = sqrt( d_a[i+1] - qwtSqr(d_b[i]));
00351         sum += qwtSqr(d_c[i]);
00352     }
00353     d_b[imax-1] = (d_b[imax-1] - d_c[imax-2] * d_b[imax-2]) / d_a[imax-1];
00354     d_a[imax] = sqrt(d_a[imax] - qwtSqr(d_b[imax-1]) - sum);
00355     
00356 
00357     // forward elimination
00358     s[0] = d[0] / d_a[0];
00359     sum = 0;
00360     for(i=1;i<imax;i++)
00361     {
00362         s[i] = (d[i] - d_b[i-1] * s[i-1]) / d_a[i];
00363         sum += d_c[i-1] * s[i-1];
00364     }
00365     s[imax] = (d[imax] - d_b[imax-1]*s[imax-1] - sum) / d_a[imax];
00366     
00367     
00368     // backward elimination
00369     s[imax] = - s[imax] / d_a[imax];
00370     s[imax-1] = -(s[imax-1] + d_b[imax-1] * s[imax]) / d_a[imax-1];
00371     for (i= imax - 2; i >= 0; i--)
00372        s[i] = - (s[i] + d_b[i] * s[i+1] + d_c[i] * s[imax]) / d_a[i];
00373 
00374     //
00375     // Finally, determine the spline coefficients
00376     //
00377     s[d_size-1] = s[0];
00378     for (i=0;i<d_size-1;i++)
00379     {
00380         d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00381         d_b[i] = 0.5 * s[i];
00382         d_c[i] = ( d_y[i+1] - d_y[i] ) 
00383             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00384     }
00385 
00386     delete[] d;
00387     delete[] s;
00388     delete[] h;
00389 
00390     return 0;
00391 }
00392 
00393 
00395 void QwtSpline::cleanup()
00396 {
00397     if (d_a) delete[] d_a;
00398     if (d_b) delete[] d_b;
00399     if (d_c) delete[] d_c;
00400     if (d_xbuffer) delete[] d_xbuffer;
00401     if (d_ybuffer) delete[] d_ybuffer;
00402     d_a = d_b = d_c = 0;
00403     d_xbuffer = d_ybuffer = d_x = d_y = 0;
00404     d_size = 0;
00405 }

Generated on Fri Nov 7 14:11:46 2003 for Qwt Developer's Guide by doxygen 1.3.2