00001
00006
00007 #include "ratint.h"
00008 #include "math.h"
00009 #include "../../Exceptions/error.h"
00010
00011 #define TINY 1.0e-25
00012
00016 void ratint(double *xa, double *ya, int n, double x, double *y, double *dy)
00017 {
00018 int m,i,ns=1;
00019 double w,t,hh,h,dd,*c,*d;
00020
00021 c = new double[n];
00022 d = new double[n];
00023 hh=fabs(x-xa[0]);
00024 for (i=0;i<=n-1;i++) {
00025 h=fabs(x-xa[i]);
00026 if (h == 0.0) {
00027 *y=ya[i];
00028 *dy=0.0;
00029 return;
00030 } else if (h < hh) {
00031 ns=i;
00032 hh=h;
00033 }
00034 c[i]=ya[i];
00035 d[i]=ya[i]+TINY;
00036 }
00037 *y=ya[ns--];
00038 for (m=0;m<n-1;m++) {
00039 for (i=0;i<n-m;i++) {
00040 w=c[i+1]-d[i];
00041 h=xa[i+m]-x;
00042 t=(xa[i]-x)*d[i]/h;
00043 dd=t-c[i+1];
00044 if (dd == 0.0) {
00045
00046
00047 throw error_msg("RATINT_ERROR", "Error in routine ratint");
00048 }
00049 dd=w/dd;
00050 d[i]=c[i+1]*dd;
00051 c[i]=t*dd;
00052 }
00053 *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
00054 }
00055 return;
00056 }