VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
ratint.cpp
Go to the documentation of this file.
1 /**
2  * Some other interpolation
3  * \file ratint.cpp
4  */
5 
6 //#include <iostream>
7 #include "ratint.h"
8 #include "math.h"
9 #include "../../Exceptions/error.h"
10 
11 #define TINY 1.0e-25
12 
13 /**
14  * Some other interpolation
15  */
16 void ratint(double *xa, double *ya, int n, double x, double *y, double *dy)
17 {
18  int m,i,ns=1;
19  double w,t,hh,h,dd,*c,*d;
20 
21  c = new double[n];//vector(1,n);
22  d = new double[n];//vector(1,n);
23  hh=fabs(x-xa[0]);
24  for (i=0;i<=n-1;i++) {
25  h=fabs(x-xa[i]);
26  if (h == 0.0) {
27  *y=ya[i];
28  *dy=0.0;
29  return;
30  } else if (h < hh) {
31  ns=i;
32  hh=h;
33  }
34  c[i]=ya[i];
35  d[i]=ya[i]+TINY;
36  }
37  *y=ya[ns--];
38  for (m=0;m<n-1;m++) {
39  for (i=0;i<n-m;i++) {
40  w=c[i+1]-d[i];
41  h=xa[i+m]-x;
42  t=(xa[i]-x)*d[i]/h;
43  dd=t-c[i+1];
44  if (dd == 0.0) {
45  //printf("Error in routine ratint");
46  //exit(0);
47  throw error_msg("RATINT_ERROR", "Error in routine ratint");
48  }
49  dd=w/dd;
50  d[i]=c[i+1]*dd;
51  c[i]=t*dd;
52  }
53  *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
54  }
55  return;
56 }