CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: DefiniteIntegral.cc,v 1.4.4.1 2003/11/21 16:48:14 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 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 } 00059 else { 00060 for (it=1,j=1;j<n-1;j++) it <<=1; 00061 double tnm=it; 00062 double del = (b-a)/tnm; 00063 double x=a+0.5*del; 00064 double sum; 00065 for (sum=0.0,j=1;j<=it;j++,x+=del) sum +=function(x); 00066 _sTrap = 0.5*(_sTrap+(b-a)*sum/tnm); 00067 } 00068 return _sTrap; 00069 } 00070 00071 00072 void DefiniteIntegral::_polint(double *xArray, double *yArray, double x, double & y, double & deltay) const { 00073 00074 double dif = fabs(x-xArray[1]),dift; 00075 double c[_KP],d[_KP]; 00076 int ns=1; 00077 for (int i=1;i<=_K;i++) { 00078 dift=fabs(x-xArray[i]); 00079 if (dift<dif) { 00080 ns=i; 00081 dif=dift; 00082 } 00083 c[i]=d[i]=yArray[i]; 00084 } 00085 y = yArray[ns--]; 00086 for (int m=1;m<_K;m++) { 00087 for (int i=1;i<=_K-m;i++) { 00088 double ho = xArray[i]-x; 00089 double hp= xArray[i+m]-x; 00090 double w=c[i+1]-d[i]; 00091 double den=ho-hp; 00092 if (den==0) 00093 std::cerr 00094 << "Error in polynomial extrapolation" 00095 << std::endl; 00096 den=w/den; 00097 d[i]=hp*den; 00098 c[i]=ho*den; 00099 } 00100 deltay = 2*ns < (_K-m) ? c[ns+1]: d[ns--]; 00101 y += deltay; 00102 } 00103 } 00104 } // namespace Genfun