00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "romberg.h"
00014
00015 namespace NR {
00016
00017 using namespace std;
00018
00019 void polint(const vector<double>& xa, const vector<double>& ya, const double x, double &y, double &dy) {
00020 int i, m, ns=0;
00021 double den, dif, dift, ho, hp, w;
00022
00023 int n = xa.size();
00024 vector<double> c(n), d(n);
00025 dif = fabs(x - xa[0]);
00026 for (i = 0; i < n; ++ i) {
00027 if ((dift = fabs(x - xa[i])) < dif) {
00028 ns = i;
00029 dif = dift;
00030 }
00031 c[i] = ya[i];
00032 d[i] = ya[i];
00033 }
00034 y = ya[ns--];
00035 for (m = 1; m < n; ++ m) {
00036 for (i = 0; i < n-m; i ++) {
00037 ho = xa[i] - x;
00038 hp = xa[i+m] - x;
00039 w = c[i+1] - d[i];
00040 den = ho - hp;
00041 assert(den != 0.0);
00042 den = w/den;
00043 d[i] = hp*den;
00044 c[i] = ho*den;
00045 }
00046 y += (dy = (2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
00047 }
00048 }
00049
00050 }