00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "qwt_spline.h"
00011 #include "qwt_math.h"
00012
00013 class QwtSpline::PrivateData
00014 {
00015 public:
00016 PrivateData():
00017 size(0),
00018 buffered(false)
00019 {
00020 a = b = c = NULL;
00021 xbuffer = ybuffer = x = y = NULL;
00022 }
00023
00024
00025 double *a;
00026 double *b;
00027 double *c;
00028 double *d;
00029
00030
00031 double *x;
00032 double *y;
00033 double *xbuffer;
00034 double *ybuffer;
00035 int size;
00036
00037
00038 bool buffered;
00039 };
00040
00042 QwtSpline::QwtSpline()
00043 {
00044 d_data = new PrivateData;
00045 }
00046
00048 QwtSpline::~QwtSpline()
00049 {
00050 cleanup();
00051 delete d_data;
00052 }
00053
00054
00077 void QwtSpline::copyValues(bool tf)
00078 {
00079 cleanup();
00080 d_data->buffered = tf;
00081 }
00082
00087 double QwtSpline::value(double x) const
00088 {
00089 if (!d_data->a)
00090 return 0.0;
00091
00092 const int i = lookup(x);
00093
00094 const double delta = x - d_data->x[i];
00095 return( ( ( ( d_data->a[i] * delta) + d_data->b[i] )
00096 * delta + d_data->c[i] ) * delta + d_data->y[i] );
00097 }
00098
00100 int QwtSpline::lookup(double x) const
00101 {
00102 int i1, i2, i3;
00103
00104 if (x <= d_data->x[0])
00105 i1 = 0;
00106 else if (x >= d_data->x[d_data->size - 2])
00107 i1 = d_data->size -2;
00108 else
00109 {
00110 i1 = 0;
00111 i2 = d_data->size -2;
00112 i3 = 0;
00113
00114 while ( i2 - i1 > 1 )
00115 {
00116 i3 = i1 + ((i2 - i1) >> 1);
00117
00118 if (d_data->x[i3] > x)
00119 i2 = i3;
00120 else
00121 i1 = i3;
00122
00123 }
00124 }
00125 return i1;
00126 }
00127
00128
00148 bool QwtSpline::recalc(double *x, double *y, int n, bool periodic)
00149 {
00150 cleanup();
00151
00152 if (n <= 2)
00153 return false;
00154
00155 d_data->size = n;
00156
00157 if (d_data->buffered)
00158 {
00159 d_data->xbuffer = new double[n];
00160 d_data->ybuffer = new double[n];
00161
00162 for (int i = 0; i < n; i++)
00163 {
00164 d_data->xbuffer[i] = x[i];
00165 d_data->ybuffer[i] = y[i];
00166 }
00167 d_data->x = d_data->xbuffer;
00168 d_data->y = d_data->ybuffer;
00169 }
00170 else
00171 {
00172 d_data->x = x;
00173 d_data->y = y;
00174 }
00175
00176 d_data->a = new double[n-1];
00177 d_data->b = new double[n-1];
00178 d_data->c = new double[n-1];
00179
00180 bool ok;
00181 if(periodic)
00182 ok = buildPerSpline();
00183 else
00184 ok = buildNatSpline();
00185
00186 if (!ok)
00187 cleanup();
00188
00189 return ok;
00190 }
00191
00196 bool QwtSpline::buildNatSpline()
00197 {
00198 int i;
00199
00200 double *d = new double[d_data->size-1];
00201 double *h = new double[d_data->size-1];
00202 double *s = new double[d_data->size];
00203
00204
00205
00206
00207 for (i = 0; i < d_data->size - 1; i++)
00208 {
00209 h[i] = d_data->x[i+1] - d_data->x[i];
00210 if (h[i] <= 0)
00211 {
00212 delete[] h;
00213 delete[] s;
00214 delete[] d;
00215 return false;
00216 }
00217 }
00218
00219 double dy1 = (d_data->y[1] - d_data->y[0]) / h[0];
00220 for (i = 1; i < d_data->size - 1; i++)
00221 {
00222 d_data->b[i] = d_data->c[i] = h[i];
00223 d_data->a[i] = 2.0 * (h[i-1] + h[i]);
00224
00225 const double dy2 = (d_data->y[i+1] - d_data->y[i]) / h[i];
00226 d[i] = 6.0 * ( dy1 - dy2);
00227 dy1 = dy2;
00228 }
00229
00230
00231
00232
00233
00234
00235 for(i = 1; i < d_data->size - 2;i++)
00236 {
00237 d_data->c[i] /= d_data->a[i];
00238 d_data->a[i+1] -= d_data->b[i] * d_data->c[i];
00239 }
00240
00241
00242 s[1] = d[1];
00243 for(i=2;i<d_data->size - 1;i++)
00244 s[i] = d[i] - d_data->c[i-1] * s[i-1];
00245
00246
00247 s[d_data->size - 2] = - s[d_data->size - 2] / d_data->a[d_data->size - 2];
00248 for (i= d_data->size -3; i > 0; i--)
00249 s[i] = - (s[i] + d_data->b[i] * s[i+1]) / d_data->a[i];
00250
00251
00252
00253
00254 s[d_data->size - 1] = s[0] = 0.0;
00255 for (i = 0; i < d_data->size - 1; i++)
00256 {
00257 d_data->a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00258 d_data->b[i] = 0.5 * s[i];
00259 d_data->c[i] = ( d_data->y[i+1] - d_data->y[i] )
00260 / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00261 }
00262
00263 delete[] d;
00264 delete[] s;
00265 delete[] h;
00266
00267 return true;
00268 }
00269
00274 bool QwtSpline::buildPerSpline()
00275 {
00276 int i;
00277
00278 double *d = new double[d_data->size-1];
00279 double *h = new double[d_data->size-1];
00280 double *s = new double[d_data->size];
00281
00282
00283
00284
00285
00286 for (i=0; i<d_data->size - 1; i++)
00287 {
00288 h[i] = d_data->x[i+1] - d_data->x[i];
00289 if (h[i] <= 0.0)
00290 {
00291 delete[] h;
00292 delete[] s;
00293 delete[] d;
00294 return false;
00295 }
00296 }
00297
00298 const int imax = d_data->size - 2;
00299 double htmp = h[imax];
00300 double dy1 = (d_data->y[0] - d_data->y[imax]) / htmp;
00301 for (i=0; i <= imax; i++)
00302 {
00303 d_data->b[i] = d_data->c[i] = h[i];
00304 d_data->a[i] = 2.0 * (htmp + h[i]);
00305 const double dy2 = (d_data->y[i+1] - d_data->y[i]) / h[i];
00306 d[i] = 6.0 * ( dy1 - dy2);
00307 dy1 = dy2;
00308 htmp = h[i];
00309 }
00310
00311
00312
00313
00314
00315
00316 d_data->a[0] = sqrt(d_data->a[0]);
00317 d_data->c[0] = h[imax] / d_data->a[0];
00318 double sum = 0;
00319
00320 for(i=0;i<imax-1;i++)
00321 {
00322 d_data->b[i] /= d_data->a[i];
00323 if (i > 0)
00324 d_data->c[i] = - d_data->c[i-1] * d_data->b[i-1] / d_data->a[i];
00325 d_data->a[i+1] = sqrt( d_data->a[i+1] - qwtSqr(d_data->b[i]));
00326 sum += qwtSqr(d_data->c[i]);
00327 }
00328 d_data->b[imax-1] = (d_data->b[imax-1] - d_data->c[imax-2] * d_data->b[imax-2]) / d_data->a[imax-1];
00329 d_data->a[imax] = sqrt(d_data->a[imax] - qwtSqr(d_data->b[imax-1]) - sum);
00330
00331
00332
00333 s[0] = d[0] / d_data->a[0];
00334 sum = 0;
00335 for(i=1;i<imax;i++)
00336 {
00337 s[i] = (d[i] - d_data->b[i-1] * s[i-1]) / d_data->a[i];
00338 sum += d_data->c[i-1] * s[i-1];
00339 }
00340 s[imax] = (d[imax] - d_data->b[imax-1]*s[imax-1] - sum) / d_data->a[imax];
00341
00342
00343
00344 s[imax] = - s[imax] / d_data->a[imax];
00345 s[imax-1] = -(s[imax-1] + d_data->b[imax-1] * s[imax]) / d_data->a[imax-1];
00346 for (i= imax - 2; i >= 0; i--)
00347 s[i] = - (s[i] + d_data->b[i] * s[i+1] + d_data->c[i] * s[imax]) / d_data->a[i];
00348
00349
00350
00351
00352 s[d_data->size-1] = s[0];
00353 for (i=0;i<d_data->size-1;i++)
00354 {
00355 d_data->a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00356 d_data->b[i] = 0.5 * s[i];
00357 d_data->c[i] = ( d_data->y[i+1] - d_data->y[i] )
00358 / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00359 }
00360
00361 delete[] d;
00362 delete[] s;
00363 delete[] h;
00364
00365 return true;
00366 }
00367
00368
00370 void QwtSpline::cleanup()
00371 {
00372 delete[] d_data->a;
00373 delete[] d_data->b;
00374 delete[] d_data->c;
00375 delete[] d_data->xbuffer;
00376 delete[] d_data->ybuffer;
00377
00378 d_data->a = d_data->b = d_data->c = NULL;
00379 d_data->xbuffer = d_data->ybuffer = d_data->x = d_data->y = NULL;
00380 d_data->size = 0;
00381 }