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
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
00214
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
00240
00241
00242
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
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
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
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
00309
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
00338
00339
00340
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
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
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
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 }