/* Slightly modified versions of the "polint" routine from * Press, William H., Brian P. Flannery, Saul A. Teukolsky and * William T. Vetterling, 1986, "Numerical Recipes: The Art of * Scientific Computing" (Fortran), Cambrigde University Press, * pp. 80-82. */ #include #include #include void polint( double *xa, double *ya, int n, double x, double *y, double *dy ); void polint( double *xa, double *ya, int n, double x, double *y, double *dy ) { double *c = NULL; double *d = NULL; double den; double dif; double dift; double ho; double hp; double w; int i; int m; int ns; if( (c = (double *)malloc( n * sizeof( double ) )) == NULL || (d = (double *)malloc( n * sizeof( double ) )) == NULL ) { fprintf( stderr, "polint error: allocating workspace\n" ); fprintf( stderr, "polint error: setting y = 0 and dy = 1e9\n" ); *y = 0.0; *dy = 1.e9; if( c != NULL ) free( c ); if( d != NULL ) free( d ); return; } ns = 0; dif = fabs(x-xa[0]); for( i = 0; i < n; ++i ) { dift = fabs( x-xa[i] ); if( dift < dif ) { ns = i; dif = dift; } c[i] = ya[i]; d[i] = ya[i]; } *y = ya[ns]; ns = ns-1; for( m = 0; m < n-1; ++m ) { for( i = 0; i < n-m-1; ++i ) { ho = xa[i]-x; hp = xa[i+m+1]-x; w = c[i+1]-d[i]; den = ho-hp; if( den == 0 ) { fprintf( stderr, "polint error: den = 0\n" ); fprintf( stderr, "polint error: setting y = 0 and dy = 1e9\n" ); *y = 0.0; *dy = 1.e9; if( c != NULL ) free( c ); if( d != NULL ) free( d ); return; } den = w/den; d[i] = hp*den; c[i] = ho*den; } if( 2*(ns+1) < n-m-1 ) { *dy = c[ns+1]; } else { *dy = d[ns]; ns = ns-1; } *y = (*y)+(*dy); } if( c != NULL ) free( c ); if( d != NULL ) free( d ); return; }