//////////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2008 The Regents of the University of California // // This file is part of Qbox // // Qbox is distributed under the terms of the GNU General Public License // as published by the Free Software Foundation, either version 2 of // the License, or (at your option) any later version. // See the file COPYING in the root directory of this distribution // or . // /////////////////////////////////////////////////////////////////////////////// // // spline.C // /////////////////////////////////////////////////////////////////////////////// // \$Id: spline.C,v 1.5 2008-09-08 15:56:20 fgygi Exp \$ #include "spline.h" #include void tridsolve(int n, double* d, double* e, double* f, double* x) { // solve the tridiagonal system Ax=b // d[i] = a(i,i) // e[i] = a(i,i+1) (superdiagonal of A, e[n-1] not defined) // f[i] = a(i,i-1) (subdiagonal of A, f[0] not defined) // x[i] = right-hand side b as input // x[i] = solution as output for ( int i = 1; i < n; i++ ) { f[i] /= d[i-1]; d[i] -= f[i]*e[i-1]; } for ( int i = 1; i < n; i++ ) x[i] -= f[i]*x[i-1]; x[n-1] /= d[n-1]; for ( int i = n-2; i >= 0; i-- ) x[i] = (x[i]-e[i]*x[i+1])/d[i]; } void spline(int n, double *x, double *y, double yp_left, double yp_right, int bcnat_left, int bcnat_right, double *y2) { const double third = 1.0/3.0; const double sixth = 1.0/6.0; double *d = new double[n]; double *e = new double[n]; double *f = new double[n]; if ( bcnat_left == 0 ) { // use derivative yp_left at x[0] const double h = x[1]-x[0]; assert(h>0.0); d[0] = third*h; e[0] = sixth*h; f[0] = 0.0; y2[0] = (y[1]-y[0])/h - yp_left; } else { // use natural spline at x[0] d[0] = 1.0; e[0] = 0.0; f[0] = 0.0; y2[0] = 0.0; } if ( bcnat_right == 0 ) { // use derivative yp_right at x[n-1] const double h = x[n-1]-x[n-2]; assert(h>0.0); d[n-1] = third*h; e[n-1] = 0.0; f[n-1] = sixth*h; y2[n-1] = yp_right - (y[n-1]-y[n-2])/h; } else { // use natural spline at x[n-1] d[n-1] = 1.0; e[n-1] = 0.0; f[n-1] = 0.0; y2[n-1] = 0.0; } // tridiagonal matrix for ( int i = 1; i < n-1; i++ ) { const double hp = x[i+1]-x[i]; const double hm = x[i]-x[i-1]; assert(hp>0.0); assert(hm>0.0); d[i] = third * (hp+hm); e[i] = sixth * hp; f[i] = sixth * hm; y2[i] = (y[i+1]-y[i])/hp - (y[i]-y[i-1])/hm; } tridsolve(n,d,e,f,y2); delete [] d; delete [] e; delete [] f; } void splint (int n, double *xa, double *ya, double *y2a, double x, double *y) { int k; double a,b,h; int kl = 0; int kh = n-1; while ( kh - kl > 1 ) { k = ( kh + kl ) / 2; if ( xa[k] > x ) kh = k; else kl = k; } h = xa[kh] - xa[kl]; assert ( h > 0.0 ); a = ( xa[kh] - x ) / h; b = ( x - xa[kl] ) / h; *y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) * ( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] ); } void splintd (int n, double *xa, double *ya, double *y2a, double x, double *y, double *dy) { int k; double a,b,h; int kl = 0; int kh = n-1; while ( kh - kl > 1 ) { k = ( kh + kl ) / 2; if ( xa[k] > x ) kh = k; else kl = k; } h = xa[kh] - xa[kl]; assert ( h > 0.0 ); a = ( xa[kh] - x ) / h; b = ( x - xa[kl] ) / h; *y = a * ya[kl] + b * ya[kh] + h * h * (1.0/6.0) * ( (a*a*a-a) * y2a[kl] + (b*b*b-b) * y2a[kh] ); *dy = ( ya[kh] - ya[kl] ) / h + h * ( ( (1.0/6.0) - 0.5 * a * a ) * y2a[kl] + ( 0.5 * b * b - (1.0/6.0) ) * y2a[kh] ); }