00001
00002
00003
00004
00005
00006
00007
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
00244
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
00270
00271
00272
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
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
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
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
00339
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
00368
00369
00370
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
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
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
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 }