00001
00006 #include "../../Logging/Output.h"
00007 #include "math.h"
00008 #include <malloc.h>
00009 #include <stdlib.h>
00010
00014 void polint( double *xa, double *ya, int n, double x, double *y, double *dy )
00015 {
00016 double *c = NULL;
00017 double *d = NULL;
00018 double den;
00019 double dif;
00020 double dift;
00021 double ho;
00022 double hp;
00023 double w;
00024
00025 int i;
00026 int m;
00027 int ns;
00028
00029 if( (c = (double *)malloc( n * sizeof( double ) )) == NULL || (d = (double *)malloc( n * sizeof( double ) )) == NULL ) {
00030
00031 Output::echo(0, "polint error: allocating workspace\n" );
00032
00033 Output::echo(0, "polint error: setting y = 0 and dy = 1e9\n" );
00034 *y = 0.0;
00035 *dy = 1.e9;
00036
00037 if( c != NULL )
00038 free( c );
00039 if( d != NULL )
00040 free( d );
00041 return;
00042 }
00043
00044 ns = 0;
00045 dif = fabs(x-xa[0]);
00046 for( i = 0; i < n; ++i ) {
00047 dift = fabs( x-xa[i] );
00048 if( dift < dif ) {
00049 ns = i;
00050 dif = dift;
00051 }
00052 c[i] = ya[i];
00053 d[i] = ya[i];
00054 }
00055 *y = ya[ns];
00056 ns = ns-1;
00057
00058 for( m = 0; m < n-1; ++m ) {
00059 for( i = 0; i < n-m-1; ++i ) {
00060 ho = xa[i]-x;
00061 hp = xa[i+m+1]-x;
00062 w = c[i+1]-d[i];
00063 den = ho-hp;
00064 if( fabs(den) < 1.e-25 ) {
00065
00066 Output::echo(0, "polint error: den = 0\n" );
00067
00068 Output::echo(0, "polint error: setting y = 0 and dy = 1e9\n" );
00069 exit(0);
00070 *y = 0.0;
00071 *dy = 1.e9;
00072
00073 if( c != NULL )
00074 free( c );
00075 if( d != NULL )
00076 free( d );
00077 return;
00078 }
00079 den = w/den;
00080 d[i] = hp*den;
00081 c[i] = ho*den;
00082 }
00083 if( 2*(ns+1) < n-m-1 ) {
00084 *dy = c[ns+1];
00085 } else {
00086 *dy = d[ns];
00087 ns = ns-1;
00088
00089 }
00090 *y = (*y)+(*dy);
00091 }
00092
00093
00094 if( c != NULL )
00095 free( c );
00096 if( d != NULL )
00097 free( d );
00098 return;
00099 }