CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: DefiniteIntegral.cc,v 1.6 2010/06/16 18:22:01 garren Exp $ 00003 00004 #include <cmath> 00005 00006 #include "CLHEP/GenericFunctions/DefiniteIntegral.hh" 00007 #include "CLHEP/GenericFunctions/AbsFunction.hh" 00008 namespace Genfun { 00009 00010 const int DefiniteIntegral::_K =5; 00011 const int DefiniteIntegral::_KP=_K+1; 00012 DefiniteIntegral::DefiniteIntegral(double a, double b): 00013 _a(a), _b(b) 00014 {} 00015 00016 DefiniteIntegral::~DefiniteIntegral() { 00017 } 00018 00019 00020 00021 #define FANCY 00022 double DefiniteIntegral::operator [] (const AbsFunction & function) const { 00023 _nFunctionCalls=0; 00024 const int MAXITERATIONS=40; // Maximum number of iterations 00025 const double EPS=1.0E-6; // Precision 00026 #ifdef FANCY 00027 00028 double s[MAXITERATIONS+2],h[MAXITERATIONS+2]; 00029 h[1]=1.0; 00030 for (int j=1;j<=MAXITERATIONS;j++) { 00031 s[j]=_trapzd(function, _a, _b, j); 00032 if (j>=_K) { 00033 double ss, dss; 00034 _polint(h+j-_K,s+j-_K,0.0,ss, dss); 00035 if (fabs(dss) <= EPS*fabs(ss)) return ss; 00036 } 00037 s[j+1]=s[j]; 00038 h[j+1]=0.25*h[j]; 00039 } 00040 #else 00041 double olds = -1E-30; 00042 for (int j=1;j<MAXITERATIONS;j++) { 00043 double s = _trapzd(function, _a,_b,j); 00044 if (fabs(s-olds)<=EPS*fabs(olds)) return s; 00045 olds=s; 00046 } 00047 #endif 00048 std::cerr 00049 << "DefiniteIntegral: too many steps. No convergence" 00050 << std::endl; 00051 return 0.0; // Never get here. 00052 } 00053 00054 double DefiniteIntegral::_trapzd(const AbsFunction & function, double a, double b, int n) const { 00055 int it, j; 00056 if (n==1) { 00057 _sTrap = 0.5*(b-a)*(function(a)+function(b)); 00058 _nFunctionCalls+=2; 00059 } 00060 else { 00061 for (it=1,j=1;j<n-1;j++) it <<=1; 00062 double tnm=it; 00063 double del = (b-a)/tnm; 00064 double x=a+0.5*del; 00065 double sum; 00066 for (sum=0.0,j=1;j<=it;j++,x+=del) { 00067 sum +=function(x); 00068 _nFunctionCalls++; 00069 } 00070 _sTrap = 0.5*(_sTrap+(b-a)*sum/tnm); 00071 } 00072 return _sTrap; 00073 } 00074 00075 00076 void DefiniteIntegral::_polint(double *xArray, double *yArray, double x, double & y, double & deltay) const { 00077 00078 double dif = fabs(x-xArray[1]),dift; 00079 double c[_KP],d[_KP]; 00080 int ns=1; 00081 for (int i=1;i<=_K;i++) { 00082 dift=fabs(x-xArray[i]); 00083 if (dift<dif) { 00084 ns=i; 00085 dif=dift; 00086 } 00087 c[i]=d[i]=yArray[i]; 00088 } 00089 y = yArray[ns--]; 00090 for (int m=1;m<_K;m++) { 00091 for (int i=1;i<=_K-m;i++) { 00092 double ho = xArray[i]-x; 00093 double hp= xArray[i+m]-x; 00094 double w=c[i+1]-d[i]; 00095 double den=ho-hp; 00096 if (den==0) 00097 std::cerr 00098 << "Error in polynomial extrapolation" 00099 << std::endl; 00100 den=w/den; 00101 d[i]=hp*den; 00102 c[i]=ho*den; 00103 } 00104 deltay = 2*ns < (_K-m) ? c[ns+1]: d[ns--]; 00105 y += deltay; 00106 } 00107 } 00108 00109 unsigned int DefiniteIntegral::numFunctionCalls() const { 00110 return _nFunctionCalls; 00111 } 00112 } // namespace Genfun