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 }