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

qwt_spline.cpp

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 
00206 int QwtSpline::recalc(const QwtArray<double> &x, const QwtArray<double> &y,
00207     int periodic)
00208 {
00209     int n = QMIN(x.size(), y.size());
00210     d_buffered = TRUE;
00211 
00212     return recalc(x.data(), y.data(), n, periodic);
00213 }
00214 
00224 int QwtSpline::buildNatSpline()
00225 {
00226     int i;
00227     double dy1, dy2;
00228     
00229     double *d = new double[d_size-1];
00230     double *h = new double[d_size-1];
00231     double *s = new double[d_size];
00232 
00233     if ( (!d) || (!h) || (!s) )
00234     {
00235         cleanup();
00236         if (h) delete[] h;
00237         if (s) delete[] s;
00238         if (d) delete[] d;
00239         return Qwt::ErrNoMem;
00240     }
00241 
00242     //
00243     //  set up tridiagonal equation system; use coefficient
00244     //  vectors as temporary buffers
00245     for (i=0; i<d_size - 1; i++) 
00246     {
00247         h[i] = d_x[i+1] - d_x[i];
00248         if (h[i] <= 0)
00249         {
00250             delete[] h;
00251             delete[] s;
00252             delete[] d;
00253             return Qwt::ErrMono;
00254         }
00255     }
00256     
00257     dy1 = (d_y[1] - d_y[0]) / h[0];
00258     for (i = 1; i < d_size - 1; i++)
00259     {
00260         d_b[i] = d_c[i] = h[i];
00261         d_a[i] = 2.0 * (h[i-1] + h[i]);
00262 
00263         dy2 = (d_y[i+1] - d_y[i]) / h[i];
00264         d[i] = 6.0 * ( dy1 - dy2);
00265         dy1 = dy2;
00266     }
00267 
00268     //
00269     // solve it
00270     //
00271     
00272     // L-U Factorization
00273     for(i = 1; i < d_size - 2;i++)
00274     {
00275         d_c[i] /= d_a[i];
00276         d_a[i+1] -= d_b[i] * d_c[i]; 
00277     }
00278 
00279     // forward elimination
00280     s[1] = d[1];
00281     for(i=2;i<d_size - 1;i++)
00282        s[i] = d[i] - d_c[i-1] * s[i-1];
00283     
00284     // backward elimination
00285     s[d_size - 2] = - s[d_size - 2] / d_a[d_size - 2];
00286     for (i= d_size -3; i > 0; i--)
00287        s[i] = - (s[i] + d_b[i] * s[i+1]) / d_a[i];
00288 
00289     //
00290     // Finally, determine the spline coefficients
00291     //
00292     s[d_size - 1] = s[0] = 0.0;
00293     for (i = 0; i < d_size - 1; i++)
00294     {
00295         d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00296         d_b[i] = 0.5 * s[i];
00297         d_c[i] = ( d_y[i+1] - d_y[i] ) 
00298             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00299     }
00300 
00301     delete[] d;
00302     delete[] s;
00303     delete[] h;
00304 
00305     return 0;
00306     
00307 }
00308 
00318 int QwtSpline::buildPerSpline()
00319 {
00320     int i,imax;
00321     double sum;
00322     double dy1, dy2,htmp;
00323     
00324     double *d = new double[d_size-1];
00325     double *h = new double[d_size-1];
00326     double *s = new double[d_size];
00327     
00328     if ( (!d) || (!h) || (!s) )
00329     {
00330         cleanup();
00331         if (h) delete[] h;
00332         if (s) delete[] s;
00333         if (d) delete[] d;
00334         return Qwt::ErrNoMem;
00335     }
00336 
00337     //
00338     //  setup equation system; use coefficient
00339     //  vectors as temporary buffers
00340     //
00341     for (i=0; i<d_size - 1; i++)
00342     {
00343         h[i] = d_x[i+1] - d_x[i];
00344         if (h[i] <= 0.0)
00345         {
00346             delete[] h;
00347             delete[] s;
00348             delete[] d;
00349             return Qwt::ErrMono;
00350         }
00351     }
00352     
00353     imax = d_size - 2;
00354     htmp = h[imax];
00355     dy1 = (d_y[0] - d_y[imax]) / htmp;
00356     for (i=0; i <= imax; i++)
00357     {
00358         d_b[i] = d_c[i] = h[i];
00359         d_a[i] = 2.0 * (htmp + h[i]);
00360         dy2 = (d_y[i+1] - d_y[i]) / h[i];
00361         d[i] = 6.0 * ( dy1 - dy2);
00362         dy1 = dy2;
00363         htmp = h[i];
00364     }
00365 
00366     //
00367     // solve it
00368     //
00369     
00370     // L-U Factorization
00371     d_a[0] = sqrt(d_a[0]);
00372     d_c[0] = h[imax] / d_a[0];
00373     sum = 0;
00374 
00375     for(i=0;i<imax-1;i++)
00376     {
00377         d_b[i] /= d_a[i];
00378         if (i > 0)
00379            d_c[i] = - d_c[i-1] * d_b[i-1] / d_a[i];
00380         d_a[i+1] = sqrt( d_a[i+1] - qwtSqr(d_b[i]));
00381         sum += qwtSqr(d_c[i]);
00382     }
00383     d_b[imax-1] = (d_b[imax-1] - d_c[imax-2] * d_b[imax-2]) / d_a[imax-1];
00384     d_a[imax] = sqrt(d_a[imax] - qwtSqr(d_b[imax-1]) - sum);
00385     
00386 
00387     // forward elimination
00388     s[0] = d[0] / d_a[0];
00389     sum = 0;
00390     for(i=1;i<imax;i++)
00391     {
00392         s[i] = (d[i] - d_b[i-1] * s[i-1]) / d_a[i];
00393         sum += d_c[i-1] * s[i-1];
00394     }
00395     s[imax] = (d[imax] - d_b[imax-1]*s[imax-1] - sum) / d_a[imax];
00396     
00397     
00398     // backward elimination
00399     s[imax] = - s[imax] / d_a[imax];
00400     s[imax-1] = -(s[imax-1] + d_b[imax-1] * s[imax]) / d_a[imax-1];
00401     for (i= imax - 2; i >= 0; i--)
00402        s[i] = - (s[i] + d_b[i] * s[i+1] + d_c[i] * s[imax]) / d_a[i];
00403 
00404     //
00405     // Finally, determine the spline coefficients
00406     //
00407     s[d_size-1] = s[0];
00408     for (i=0;i<d_size-1;i++)
00409     {
00410         d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00411         d_b[i] = 0.5 * s[i];
00412         d_c[i] = ( d_y[i+1] - d_y[i] ) 
00413             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00414     }
00415 
00416     delete[] d;
00417     delete[] s;
00418     delete[] h;
00419 
00420     return 0;
00421 }
00422 
00423 
00425 void QwtSpline::cleanup()
00426 {
00427     if (d_a) delete[] d_a;
00428     if (d_b) delete[] d_b;
00429     if (d_c) delete[] d_c;
00430     if (d_xbuffer) delete[] d_xbuffer;
00431     if (d_ybuffer) delete[] d_ybuffer;
00432     d_a = d_b = d_c = 0;
00433     d_xbuffer = d_ybuffer = d_x = d_y = 0;
00434     d_size = 0;
00435 }

Generated on Sun Nov 21 11:12:44 2004 for Qwt User's Guide by doxygen 1.3.5