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

DefiniteIntegral.cc
Go to the documentation of this file.
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