10 const int DefiniteIntegral::_K =5;
11 const int DefiniteIntegral::_KP=_K+1;
24 const int MAXITERATIONS=40;
25 const double EPS=1.0E-6;
28 double s[MAXITERATIONS+2],h[MAXITERATIONS+2];
30 for (
int j=1;j<=MAXITERATIONS;j++) {
31 s[j]=_trapzd(
function, _a, _b, j);
34 _polint(h+j-_K,s+j-_K,0.0,ss, dss);
35 if (fabs(dss) <= EPS*fabs(ss))
return ss;
42 for (
int j=1;j<MAXITERATIONS;j++) {
43 double s = _trapzd(
function, _a,_b,j);
44 if (fabs(s-olds)<=EPS*fabs(olds))
return s;
49 <<
"DefiniteIntegral: too many steps. No convergence"
54 double DefiniteIntegral::_trapzd(
const AbsFunction &
function,
double a,
double b,
int n)
const {
57 _sTrap = 0.5*(b-
a)*(
function(a)+
function(
b));
61 for (it=1,j=1;j<n-1;j++) it <<=1;
63 double del = (b-
a)/tnm;
66 for (sum=0.0,j=1;j<=it;j++,x+=del) {
70 _sTrap = 0.5*(_sTrap+(b-
a)*sum/tnm);
76 void DefiniteIntegral::_polint(
double *xArray,
double *yArray,
double x,
double & y,
double & deltay)
const {
78 double dif = fabs(x-xArray[1]),dift;
81 for (
int i=1;i<=_K;i++) {
82 dift=fabs(x-xArray[i]);
90 for (
int m=1;
m<_K;
m++) {
91 for (
int i=1;i<=_K-
m;i++) {
92 double ho = xArray[i]-x;
93 double hp= xArray[i+m]-x;
98 <<
"Error in polynomial extrapolation"
104 deltay = 2*ns < (_K-m) ? c[ns+1]: d[ns--];
110 return _nFunctionCalls;