CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

DefiniteIntegral.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: DefiniteIntegral.cc,v 1.6 2010/06/16 18:22:01 garren Exp $
3 
4 #include <cmath>
5 
8 namespace Genfun {
9 
10 const int DefiniteIntegral::_K =5;
11 const int DefiniteIntegral::_KP=_K+1;
13 _a(a), _b(b)
14 {}
15 
17 }
18 
19 
20 
21 #define FANCY
22 double DefiniteIntegral::operator [] (const AbsFunction & function) const {
23  _nFunctionCalls=0;
24  const int MAXITERATIONS=40; // Maximum number of iterations
25  const double EPS=1.0E-6; // Precision
26 #ifdef FANCY
27 
28  double s[MAXITERATIONS+2],h[MAXITERATIONS+2];
29  h[1]=1.0;
30  for (int j=1;j<=MAXITERATIONS;j++) {
31  s[j]=_trapzd(function, _a, _b, j);
32  if (j>=_K) {
33  double ss, dss;
34  _polint(h+j-_K,s+j-_K,0.0,ss, dss);
35  if (fabs(dss) <= EPS*fabs(ss)) return ss;
36  }
37  s[j+1]=s[j];
38  h[j+1]=0.25*h[j];
39  }
40 #else
41  double olds = -1E-30;
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;
45  olds=s;
46  }
47 #endif
48  std::cerr
49  << "DefiniteIntegral: too many steps. No convergence"
50  << std::endl;
51  return 0.0; // Never get here.
52 }
53 
54 double DefiniteIntegral::_trapzd(const AbsFunction & function, double a, double b, int n) const {
55  int it, j;
56  if (n==1) {
57  _sTrap = 0.5*(b-a)*(function(a)+function(b));
58  _nFunctionCalls+=2;
59  }
60  else {
61  for (it=1,j=1;j<n-1;j++) it <<=1;
62  double tnm=it;
63  double del = (b-a)/tnm;
64  double x=a+0.5*del;
65  double sum;
66  for (sum=0.0,j=1;j<=it;j++,x+=del) {
67  sum +=function(x);
68  _nFunctionCalls++;
69  }
70  _sTrap = 0.5*(_sTrap+(b-a)*sum/tnm);
71  }
72  return _sTrap;
73 }
74 
75 
76 void DefiniteIntegral::_polint(double *xArray, double *yArray, double x, double & y, double & deltay) const {
77 
78  double dif = fabs(x-xArray[1]),dift;
79  double c[_KP],d[_KP];
80  int ns=1;
81  for (int i=1;i<=_K;i++) {
82  dift=fabs(x-xArray[i]);
83  if (dift<dif) {
84  ns=i;
85  dif=dift;
86  }
87  c[i]=d[i]=yArray[i];
88  }
89  y = yArray[ns--];
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;
94  double w=c[i+1]-d[i];
95  double den=ho-hp;
96  if (den==0)
97  std::cerr
98  << "Error in polynomial extrapolation"
99  << std::endl;
100  den=w/den;
101  d[i]=hp*den;
102  c[i]=ho*den;
103  }
104  deltay = 2*ns < (_K-m) ? c[ns+1]: d[ns--];
105  y += deltay;
106  }
107 }
108 
109  unsigned int DefiniteIntegral::numFunctionCalls() const {
110  return _nFunctionCalls;
111  }
112 } // namespace Genfun