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 Tue Nov 16 21:12:21 2004 for Qwt User's Guide by doxygen 1.3.8