G4Integrator< T, F > Class Template Reference

#include <G4Integrator.hh>


Public Member Functions

 G4Integrator ()
 ~G4Integrator ()
G4double Simpson (T &typeT, F f, G4double a, G4double b, G4int n)
G4double Simpson (T *ptrT, F f, G4double a, G4double b, G4int n)
G4double Simpson (G4double(*f)(G4double), G4double a, G4double b, G4int n)
G4double AdaptiveGauss (T &typeT, F f, G4double a, G4double b, G4double e)
G4double AdaptiveGauss (T *ptrT, F f, G4double a, G4double b, G4double e)
G4double AdaptiveGauss (G4double(*f)(G4double), G4double a, G4double b, G4double e)
G4double Legendre (T &typeT, F f, G4double a, G4double b, G4int n)
G4double Legendre (T *ptrT, F f, G4double a, G4double b, G4int n)
G4double Legendre (G4double(*f)(G4double), G4double a, G4double b, G4int n)
G4double Legendre10 (T &typeT, F f, G4double a, G4double b)
G4double Legendre10 (T *ptrT, F f, G4double a, G4double b)
G4double Legendre10 (G4double(*f)(G4double), G4double a, G4double b)
G4double Legendre96 (T &typeT, F f, G4double a, G4double b)
G4double Legendre96 (T *ptrT, F f, G4double a, G4double b)
G4double Legendre96 (G4double(*f)(G4double), G4double a, G4double b)
G4double Chebyshev (T &typeT, F f, G4double a, G4double b, G4int n)
G4double Chebyshev (T *ptrT, F f, G4double a, G4double b, G4int n)
G4double Chebyshev (G4double(*f)(G4double), G4double a, G4double b, G4int n)
G4double Laguerre (T &typeT, F f, G4double alpha, G4int n)
G4double Laguerre (T *ptrT, F f, G4double alpha, G4int n)
G4double Laguerre (G4double(*f)(G4double), G4double alpha, G4int n)
G4double Hermite (T &typeT, F f, G4int n)
G4double Hermite (T *ptrT, F f, G4int n)
G4double Hermite (G4double(*f)(G4double), G4int n)
G4double Jacobi (T &typeT, F f, G4double alpha, G4double beta, G4int n)
G4double Jacobi (T *ptrT, F f, G4double alpha, G4double beta, G4int n)
G4double Jacobi (G4double(*f)(G4double), G4double alpha, G4double beta, G4int n)

Protected Member Functions

G4double Gauss (T &typeT, F f, G4double a, G4double b)
G4double Gauss (T *ptrT, F f, G4double a, G4double b)
G4double Gauss (G4double(*f)(G4double), G4double a, G4double b)
void AdaptGauss (T &typeT, F f, G4double a, G4double b, G4double e, G4double &sum, G4int &n)
void AdaptGauss (T *typeT, F f, G4double a, G4double b, G4double e, G4double &sum, G4int &n)
void AdaptGauss (G4double(*f)(G4double), G4double a, G4double b, G4double e, G4double &sum, G4int &n)
G4double GammaLogarithm (G4double xx)


Detailed Description

template<class T, class F>
class G4Integrator< T, F >

Definition at line 49 of file G4Integrator.hh.


Constructor & Destructor Documentation

template<class T, class F>
G4Integrator< T, F >::G4Integrator (  )  [inline]

Definition at line 53 of file G4Integrator.hh.

00053 {;}

template<class T, class F>
G4Integrator< T, F >::~G4Integrator (  )  [inline]

Definition at line 54 of file G4Integrator.hh.

00054 {;}


Member Function Documentation

template<class T, class F>
void G4Integrator< T, F >::AdaptGauss ( G4double(*)(G4double f,
G4double  a,
G4double  b,
G4double  e,
G4double sum,
G4int n 
) [protected]

Definition at line 227 of file G4Integrator.icc.

References G4Integrator< T, F >::AdaptGauss(), G4cout, G4endl, and G4Integrator< T, F >::Gauss().

00231 {
00232    if(depth > 100)
00233    {
00234      G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl  ;
00235      G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
00236            <<G4endl ;
00237 
00238      return ;
00239    }
00240    G4double xMean = (xInitial + xFinal)/2.0 ;
00241    G4double leftHalf  = Gauss(f,xInitial,xMean) ;
00242    G4double rightHalf = Gauss(f,xMean,xFinal) ;
00243    G4double full = Gauss(f,xInitial,xFinal) ;
00244    if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
00245    {
00246       sum += full ;
00247    }
00248    else
00249    {
00250       depth++ ;
00251       AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ;
00252       AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ;
00253    }
00254 }

template<class T, class F>
void G4Integrator< T, F >::AdaptGauss ( T *  typeT,
f,
G4double  a,
G4double  b,
G4double  e,
G4double sum,
G4int n 
) [protected]

Definition at line 215 of file G4Integrator.icc.

References G4Integrator< T, F >::AdaptGauss().

00219 {
00220   AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ;
00221 }

template<class T, class F>
void G4Integrator< T, F >::AdaptGauss ( T &  typeT,
f,
G4double  a,
G4double  b,
G4double  e,
G4double sum,
G4int n 
) [protected]

Definition at line 185 of file G4Integrator.icc.

References G4cout, G4endl, and G4Integrator< T, F >::Gauss().

Referenced by G4Integrator< T, F >::AdaptGauss(), and G4Integrator< T, F >::AdaptiveGauss().

00189 {
00190    if(depth > 100)
00191    {
00192      G4cout<<"G4Integrator<T,F>::AdaptGauss: WARNING !!!"<<G4endl  ;
00193      G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps "
00194            <<G4endl ;
00195 
00196      return ;
00197    }
00198    G4double xMean = (xInitial + xFinal)/2.0 ;
00199    G4double leftHalf  = Gauss(typeT,f,xInitial,xMean) ;
00200    G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ;
00201    G4double full = Gauss(typeT,f,xInitial,xFinal) ;
00202    if(std::fabs(leftHalf+rightHalf-full) < fTolerance)
00203    {
00204       sum += full ;
00205    }
00206    else
00207    {
00208       depth++ ;
00209       AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ;
00210       AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ;
00211    }
00212 }

template<class T, class F>
G4double G4Integrator< T, F >::AdaptiveGauss ( G4double(*)(G4double f,
G4double  a,
G4double  b,
G4double  e 
)

Definition at line 289 of file G4Integrator.icc.

References G4Integrator< T, F >::AdaptGauss().

00291 {
00292    G4int depth = 0 ;
00293    G4double sum = 0.0 ;
00294    AdaptGauss(f,xInitial,xFinal,e,sum,depth) ;
00295    return sum ;
00296 }

template<class T, class F>
G4double G4Integrator< T, F >::AdaptiveGauss ( T *  ptrT,
f,
G4double  a,
G4double  b,
G4double  e 
)

Definition at line 277 of file G4Integrator.icc.

References G4Integrator< T, F >::AdaptiveGauss().

00279 {
00280   return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ;
00281 }

template<class T, class F>
G4double G4Integrator< T, F >::AdaptiveGauss ( T &  typeT,
f,
G4double  a,
G4double  b,
G4double  e 
)

Definition at line 262 of file G4Integrator.icc.

References G4Integrator< T, F >::AdaptGauss().

Referenced by G4Integrator< T, F >::AdaptiveGauss(), G4NuclNuclDiffuseElastic::TestAngleTable(), and G4DiffuseElastic::TestAngleTable().

00264 {
00265    G4int depth = 0 ;
00266    G4double sum = 0.0 ;
00267    AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ;
00268    return sum ;
00269 }

template<class T, class F>
G4double G4Integrator< T, F >::Chebyshev ( G4double(*)(G4double f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 794 of file G4Integrator.icc.

References G4INCL::Math::pi.

00796 {
00797    G4int i ;
00798    G4double xDiff, xMean, dx, integral = 0.0 ;
00799    
00800    G4int fNumber = nChebyshev  ;   // Try to reduce fNumber twice ??
00801    G4double cof = CLHEP::pi/fNumber ;
00802    G4double* fAbscissa = new G4double[fNumber] ;
00803    G4double* fWeight   = new G4double[fNumber] ;
00804    for(i=0;i<fNumber;i++)
00805    {
00806       fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
00807       fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
00808    }
00809 
00810    //
00811    // Now we ready to estimate the integral
00812    //
00813 
00814    xMean = 0.5*(a + b) ;
00815    xDiff = 0.5*(b - a) ;
00816    for(i=0;i<fNumber;i++)
00817    {
00818       dx = xDiff*fAbscissa[i] ;
00819       integral += fWeight[i]*(*f)(xMean + dx)  ;
00820    }
00821    delete[] fAbscissa;
00822    delete[] fWeight;
00823    return integral *= xDiff ;
00824 }

template<class T, class F>
G4double G4Integrator< T, F >::Chebyshev ( T *  ptrT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 783 of file G4Integrator.icc.

References G4Integrator< T, F >::Chebyshev().

00785 {
00786   return Chebyshev(*ptrT,f,a,b,n) ;
00787 } 

template<class T, class F>
G4double G4Integrator< T, F >::Chebyshev ( T &  typeT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 746 of file G4Integrator.icc.

References G4INCL::Math::pi.

Referenced by G4Integrator< T, F >::Chebyshev().

00748 {
00749    G4int i ;
00750    G4double xDiff, xMean, dx, integral = 0.0 ;
00751    
00752    G4int fNumber = nChebyshev  ;   // Try to reduce fNumber twice ??
00753    G4double cof = CLHEP::pi/fNumber ;
00754    G4double* fAbscissa = new G4double[fNumber] ;
00755    G4double* fWeight   = new G4double[fNumber] ;
00756    for(i=0;i<fNumber;i++)
00757    {
00758       fAbscissa[i] = std::cos(cof*(i + 0.5)) ;
00759       fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ;
00760    }
00761 
00762    //
00763    // Now we ready to estimate the integral
00764    //
00765 
00766    xMean = 0.5*(a + b) ;
00767    xDiff = 0.5*(b - a) ;
00768    for(i=0;i<fNumber;i++)
00769    {
00770       dx = xDiff*fAbscissa[i] ;
00771       integral += fWeight[i]*(typeT.*f)(xMean + dx)  ;
00772    }
00773    delete[] fAbscissa;
00774    delete[] fWeight;
00775    return integral *= xDiff ;
00776 }

template<class T, class F>
G4double G4Integrator< T, F >::GammaLogarithm ( G4double  xx  )  [protected]

Definition at line 1016 of file G4Integrator.icc.

Referenced by G4Integrator< T, F >::Jacobi(), and G4Integrator< T, F >::Laguerre().

01017 {
01018   static G4double cof[6] = { 76.18009172947146,     -86.50532032941677,
01019                              24.01409824083091,      -1.231739572450155,
01020                               0.1208650973866179e-2, -0.5395239384953e-5  } ;
01021   register G4int j;
01022   G4double x = xx - 1.0 ;
01023   G4double tmp = x + 5.5 ;
01024   tmp -= (x + 0.5) * std::log(tmp) ;
01025   G4double ser = 1.000000000190015 ;
01026 
01027   for ( j = 0; j <= 5; j++ )
01028   {
01029     x += 1.0 ;
01030     ser += cof[j]/x ;
01031   }
01032   return -tmp + std::log(2.5066282746310005*ser) ;
01033 }

template<class T, class F>
G4double G4Integrator< T, F >::Gauss ( G4double(*)(G4double f,
G4double  a,
G4double  b 
) [protected]

Definition at line 167 of file G4Integrator.icc.

00169 {
00170    static G4double root = 1.0/std::sqrt(3.0) ;
00171    
00172    G4double xMean = (xInitial + xFinal)/2.0 ;
00173    G4double Step  = (xFinal - xInitial)/2.0 ;
00174    G4double delta = Step*root ;
00175    G4double sum   = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ;
00176    
00177    return sum*Step ;   
00178 }

template<class T, class F>
G4double G4Integrator< T, F >::Gauss ( T *  ptrT,
f,
G4double  a,
G4double  b 
) [protected]

Definition at line 157 of file G4Integrator.icc.

References G4Integrator< T, F >::Gauss().

00158 {
00159   return Gauss(*ptrT,f,a,b) ;
00160 }

template<class T, class F>
G4double G4Integrator< T, F >::Gauss ( T &  typeT,
f,
G4double  a,
G4double  b 
) [protected]

Definition at line 138 of file G4Integrator.icc.

Referenced by G4Integrator< T, F >::AdaptGauss(), and G4Integrator< T, F >::Gauss().

00140 {
00141    static G4double root = 1.0/std::sqrt(3.0) ;
00142    
00143    G4double xMean = (xInitial + xFinal)/2.0 ;
00144    G4double Step = (xFinal - xInitial)/2.0 ;
00145    G4double delta = Step*root ;
00146    G4double sum = ((typeT.*f)(xMean + delta) + 
00147                    (typeT.*f)(xMean - delta)) ;
00148    
00149    return sum*Step ;   
00150 }

template<class T, class F>
G4double G4Integrator< T, F >::Hermite ( G4double(*)(G4double f,
G4int  n 
)

Definition at line 1145 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4INCL::Math::pi.

01146 {
01147    const G4double tolerance = 1.0e-12 ;
01148    const G4int maxNumber = 12 ;
01149    
01150    G4int i, j, k ;
01151    G4double integral = 0.0 ;
01152    G4double nwt=0., nwt1, temp1, temp2, temp3, temp ;
01153 
01154    G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ;    // 1.0/std::sqrt(std::sqrt(pi)) ??
01155 
01156    G4int fNumber = (nHermite +1)/2 ;
01157    G4double* fAbscissa = new G4double[fNumber] ;
01158    G4double* fWeight   = new G4double[fNumber] ;
01159 
01160    for(i=1;i<=fNumber;i++)
01161    {
01162       if(i == 1)
01163       {
01164          nwt = std::sqrt((G4double)(2*nHermite + 1)) - 
01165                1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
01166       }
01167       else if(i == 2)
01168       {
01169          nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ;
01170       }
01171       else if(i == 3)
01172       {
01173          nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ;
01174       }
01175       else if(i == 4)
01176       {
01177          nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ;
01178       }
01179       else 
01180       {
01181          nwt = 2.0*nwt - fAbscissa[i - 3] ;
01182       }
01183       for(k=1;k<=maxNumber;k++)
01184       {
01185          temp1 = piInMinusQ ;
01186          temp2 = 0.0 ;
01187 
01188          for(j=1;j<=nHermite;j++)
01189          {
01190             temp3 = temp2 ;
01191             temp2 = temp1 ;
01192             temp1 = nwt*std::sqrt(2.0/j)*temp2 - 
01193                     std::sqrt(((G4double)(j - 1))/j)*temp3 ;
01194          }
01195          temp = std::sqrt((G4double)2*nHermite)*temp2 ;
01196          nwt1 = nwt ;
01197          nwt = nwt1 - temp1/temp ;
01198 
01199          if(std::fabs(nwt - nwt1) <= tolerance) 
01200          {
01201             break ;
01202          }
01203       }
01204       if(k > maxNumber)
01205       {
01206          G4Exception("G4Integrator<T,F>::Hermite(...)", "Error",
01207                      FatalException, "Too many (>12) iterations.");
01208       }
01209       fAbscissa[i-1] =  nwt ;
01210       fWeight[i-1] = 2.0/(temp*temp) ;
01211    }
01212 
01213    //
01214    // Integral calculation
01215    //
01216 
01217    for(i=0;i<fNumber;i++)
01218    {
01219      integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i])   ) ;
01220    }
01221    delete[] fAbscissa;
01222    delete[] fWeight;
01223    return integral ;
01224 }

template<class T, class F>
G4double G4Integrator< T, F >::Hermite ( T *  ptrT,
f,
G4int  n 
)

Definition at line 1135 of file G4Integrator.icc.

References G4Integrator< T, F >::Hermite().

01136 {
01137   return Hermite(*ptrT,f,n) ;
01138 } 

template<class T, class F>
G4double G4Integrator< T, F >::Hermite ( T &  typeT,
f,
G4int  n 
)

Definition at line 1047 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4INCL::Math::pi.

Referenced by G4Integrator< T, F >::Hermite().

01048 {
01049    const G4double tolerance = 1.0e-12 ;
01050    const G4int maxNumber = 12 ;
01051    
01052    G4int i, j, k ;
01053    G4double integral = 0.0 ;
01054    G4double nwt=0., nwt1, temp1, temp2, temp3, temp ;
01055 
01056    G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
01057 
01058    G4int fNumber = (nHermite +1)/2 ;
01059    G4double* fAbscissa = new G4double[fNumber] ;
01060    G4double* fWeight   = new G4double[fNumber] ;
01061 
01062    for(i=1;i<=fNumber;i++)
01063    {
01064       if(i == 1)
01065       {
01066          nwt = std::sqrt((G4double)(2*nHermite + 1)) - 
01067                1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
01068       }
01069       else if(i == 2)
01070       {
01071          nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ;
01072       }
01073       else if(i == 3)
01074       {
01075          nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ;
01076       }
01077       else if(i == 4)
01078       {
01079          nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ;
01080       }
01081       else 
01082       {
01083          nwt = 2.0*nwt - fAbscissa[i - 3] ;
01084       }
01085       for(k=1;k<=maxNumber;k++)
01086       {
01087          temp1 = piInMinusQ ;
01088          temp2 = 0.0 ;
01089 
01090          for(j=1;j<=nHermite;j++)
01091          {
01092             temp3 = temp2 ;
01093             temp2 = temp1 ;
01094             temp1 = nwt*std::sqrt(2.0/j)*temp2 - 
01095                     std::sqrt(((G4double)(j - 1))/j)*temp3 ;
01096          }
01097          temp = std::sqrt((G4double)2*nHermite)*temp2 ;
01098          nwt1 = nwt ;
01099          nwt = nwt1 - temp1/temp ;
01100 
01101          if(std::fabs(nwt - nwt1) <= tolerance) 
01102          {
01103             break ;
01104          }
01105       }
01106       if(k > maxNumber)
01107       {
01108          G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error",
01109                      FatalException, "Too many (>12) iterations.");
01110       }
01111       fAbscissa[i-1] =  nwt ;
01112       fWeight[i-1] = 2.0/(temp*temp) ;
01113    }
01114 
01115    //
01116    // Integral calculation
01117    //
01118 
01119    for(i=0;i<fNumber;i++)
01120    {
01121      integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) + 
01122                               (typeT.*f)(-fAbscissa[i])   ) ;
01123    }
01124    delete[] fAbscissa;
01125    delete[] fWeight;
01126    return integral ;
01127 }

template<class T, class F>
G4double G4Integrator< T, F >::Jacobi ( G4double(*)(G4double f,
G4double  alpha,
G4double  beta,
G4int  n 
)

Definition at line 1365 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().

01367 {
01368   const G4double tolerance = 1.0e-12 ;
01369   const G4double maxNumber = 12 ;
01370   G4int i, k, j ;
01371   G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
01372   G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ;
01373 
01374   G4int     fNumber   = nJacobi ;
01375   G4double* fAbscissa = new G4double[fNumber] ;
01376   G4double* fWeight   = new G4double[fNumber] ;
01377 
01378   for (i=1;i<=nJacobi;i++)
01379   {
01380      if (i == 1)
01381      {
01382         alphaReduced = alpha/nJacobi ;
01383         betaReduced = beta/nJacobi ;
01384         root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
01385               0.767999*alphaReduced/nJacobi) ;
01386         root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
01387                 0.451998*alphaReduced*alphaReduced +
01388                 0.83001*alphaReduced*betaReduced      ;
01389         root  = 1.0-root1/root2 ;
01390      } 
01391      else if (i == 2)
01392      {
01393         root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
01394         root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
01395         root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
01396         root -= (1.0-root)*root1*root2*root3 ;
01397      } 
01398      else if (i == 3) 
01399      {
01400         root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
01401         root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
01402         root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
01403         root -= (fAbscissa[0]-root)*root1*root2*root3 ;
01404      }
01405      else if (i == nJacobi-1)
01406      {
01407         root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
01408         root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
01409         root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
01410         root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
01411      } 
01412      else if (i == nJacobi) 
01413      {
01414         root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
01415         root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
01416         root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
01417         root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
01418      } 
01419      else
01420      {
01421         root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
01422      }
01423      alphaBeta = alpha + beta ;
01424      for (k=1;k<=maxNumber;k++)
01425      {
01426         temp = 2.0 + alphaBeta ;
01427         nwt1 = (alpha-beta+temp*root)/2.0 ;
01428         nwt2 = 1.0 ;
01429         for (j=2;j<=nJacobi;j++)
01430         {
01431            nwt3 = nwt2 ;
01432            nwt2 = nwt1 ;
01433            temp = 2*j+alphaBeta ;
01434            a = 2*j*(j+alphaBeta)*(temp-2.0) ;
01435            b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
01436            c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
01437            nwt1 = (b*nwt2-c*nwt3)/a ;
01438         }
01439         nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 +
01440              2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2) /
01441              (temp*(1.0 - root*root)) ;
01442         rootTemp = root ;
01443         root = rootTemp - nwt1/nwt ;
01444         if (std::fabs(root-rootTemp) <= tolerance)
01445         {
01446            break ;
01447         }
01448      }
01449      if (k > maxNumber) 
01450      {
01451         G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error",
01452                     FatalException, "Too many (>12) iterations.");
01453      }
01454      fAbscissa[i-1] = root ;
01455      fWeight[i-1] =
01456         std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 
01457                  GammaLogarithm((G4double)(beta+nJacobi)) - 
01458                  GammaLogarithm((G4double)(nJacobi+1.0)) -
01459                  GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
01460         *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2);
01461    }
01462 
01463    //
01464    // Calculation of the integral
01465    //
01466 
01467    G4double integral = 0.0 ;
01468    for(i=0;i<fNumber;i++)
01469    {
01470       integral += fWeight[i]*(*f)(fAbscissa[i]) ;
01471    }
01472    delete[] fAbscissa;
01473    delete[] fWeight;
01474    return integral ;
01475 }

template<class T, class F>
G4double G4Integrator< T, F >::Jacobi ( T *  ptrT,
f,
G4double  alpha,
G4double  beta,
G4int  n 
)

Definition at line 1354 of file G4Integrator.icc.

References G4Integrator< T, F >::Jacobi().

01356 {
01357   return Jacobi(*ptrT,f,alpha,beta,n) ;
01358 } 

template<class T, class F>
G4double G4Integrator< T, F >::Jacobi ( T &  typeT,
f,
G4double  alpha,
G4double  beta,
G4int  n 
)

Definition at line 1237 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().

Referenced by G4Integrator< T, F >::Jacobi().

01239 {
01240   const G4double tolerance = 1.0e-12 ;
01241   const G4double maxNumber = 12 ;
01242   G4int i, k, j ;
01243   G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ;
01244   G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ;
01245 
01246   G4int     fNumber   = nJacobi ;
01247   G4double* fAbscissa = new G4double[fNumber] ;
01248   G4double* fWeight   = new G4double[fNumber] ;
01249 
01250   for (i=1;i<=nJacobi;i++)
01251   {
01252      if (i == 1)
01253      {
01254         alphaReduced = alpha/nJacobi ;
01255         betaReduced = beta/nJacobi ;
01256         root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+
01257               0.767999*alphaReduced/nJacobi) ;
01258         root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced +
01259                 0.451998*alphaReduced*alphaReduced +
01260                 0.83001*alphaReduced*betaReduced      ;
01261         root  = 1.0-root1/root2 ;
01262      } 
01263      else if (i == 2)
01264      {
01265         root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ;
01266         root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ;
01267         root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ;
01268         root -= (1.0-root)*root1*root2*root3 ;
01269      } 
01270      else if (i == 3) 
01271      {
01272         root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ;
01273         root2=1.0+0.22*(nJacobi-8.0)/nJacobi ;
01274         root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ;
01275         root -= (fAbscissa[0]-root)*root1*root2*root3 ;
01276      }
01277      else if (i == nJacobi-1)
01278      {
01279         root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ;
01280         root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ;
01281         root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ;
01282         root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ;
01283      } 
01284      else if (i == nJacobi) 
01285      {
01286         root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ;
01287         root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ;
01288         root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ;
01289         root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ;
01290      } 
01291      else
01292      {
01293         root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ;
01294      }
01295      alphaBeta = alpha + beta ;
01296      for (k=1;k<=maxNumber;k++)
01297      {
01298         temp = 2.0 + alphaBeta ;
01299         nwt1 = (alpha-beta+temp*root)/2.0 ;
01300         nwt2 = 1.0 ;
01301         for (j=2;j<=nJacobi;j++)
01302         {
01303            nwt3 = nwt2 ;
01304            nwt2 = nwt1 ;
01305            temp = 2*j+alphaBeta ;
01306            a = 2*j*(j+alphaBeta)*(temp-2.0) ;
01307             b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ;
01308            c = 2.0*(j-1+alpha)*(j-1+beta)*temp ;
01309            nwt1 = (b*nwt2-c*nwt3)/a ;
01310         }
01311         nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 +
01312               2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2)/
01313              (temp*(1.0 - root*root)) ;
01314         rootTemp = root ;
01315         root = rootTemp - nwt1/nwt ;
01316         if (std::fabs(root-rootTemp) <= tolerance)
01317         {
01318            break ;
01319         }
01320      }
01321      if (k > maxNumber) 
01322      {
01323         G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error",
01324                     FatalException, "Too many (>12) iterations.");
01325      }
01326      fAbscissa[i-1] = root ;
01327      fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 
01328                         GammaLogarithm((G4double)(beta+nJacobi)) - 
01329                         GammaLogarithm((G4double)(nJacobi+1.0)) -
01330                         GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0)))
01331                         *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2)             ;
01332    }
01333 
01334    //
01335    // Calculation of the integral
01336    //
01337 
01338    G4double integral = 0.0 ;
01339    for(i=0;i<fNumber;i++)
01340    {
01341       integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
01342    }
01343    delete[] fAbscissa;
01344    delete[] fWeight;
01345    return integral ;
01346 }

template<class T, class F>
G4double G4Integrator< T, F >::Laguerre ( G4double(*)(G4double f,
G4double  alpha,
G4int  n 
)

Definition at line 932 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().

00934 {
00935    const G4double tolerance = 1.0e-10 ;
00936    const G4int maxNumber = 12 ;
00937    G4int i, j, k ;
00938    G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ;
00939    G4double integral = 0.0 ;
00940 
00941    G4int fNumber = nLaguerre ;
00942    G4double* fAbscissa = new G4double[fNumber] ;
00943    G4double* fWeight   = new G4double[fNumber] ;
00944       
00945    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
00946    {
00947       if(i == 1)
00948       {
00949          nwt = (1.0 + alpha)*(3.0 + 0.92*alpha)
00950              / (1.0 + 2.4*fNumber + 1.8*alpha) ;
00951       }
00952       else if(i == 2)
00953       {
00954          nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
00955       }
00956       else
00957       {
00958          cofi = i - 2 ;
00959          nwt += ((1.0+2.55*cofi)/(1.9*cofi)
00960               + 1.26*cofi*alpha/(1.0+3.5*cofi))
00961               * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
00962       }
00963       for(k=1;k<=maxNumber;k++)
00964       {
00965          temp1 = 1.0 ;
00966          temp2 = 0.0 ;
00967 
00968          for(j=1;j<=fNumber;j++)
00969          {
00970             temp3 = temp2 ;
00971             temp2 = temp1 ;
00972          temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ;
00973          }
00974          temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ;
00975          nwt1 = nwt ;
00976          nwt  = nwt1 - temp1/temp ;
00977 
00978          if(std::fabs(nwt - nwt1) <= tolerance) 
00979          {
00980             break ;
00981          }
00982       }
00983       if(k > maxNumber)
00984       {
00985          G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error",
00986                      FatalException, "Too many (>12) iterations.");
00987       }
00988          
00989       fAbscissa[i-1] =  nwt ;
00990       fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 
00991                 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
00992    }
00993 
00994    //
00995    // Integral evaluation
00996    //
00997 
00998    for(i=0;i<fNumber;i++)
00999    {
01000       integral += fWeight[i]*(*f)(fAbscissa[i]) ;
01001    }
01002    delete[] fAbscissa;
01003    delete[] fWeight;
01004    return integral ;
01005 }

template<class T, class F>
G4double G4Integrator< T, F >::Laguerre ( T *  ptrT,
f,
G4double  alpha,
G4int  n 
)

Definition at line 922 of file G4Integrator.icc.

References G4Integrator< T, F >::Laguerre().

00923 {
00924   return Laguerre(*ptrT,f,alpha,nLaguerre) ;
00925 }

template<class T, class F>
G4double G4Integrator< T, F >::Laguerre ( T &  typeT,
f,
G4double  alpha,
G4int  n 
)

Definition at line 840 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().

Referenced by G4SynchrotronRadiationInMat::GetAngleK(), G4SynchrotronRadiationInMat::GetEnergyProbSR(), G4SynchrotronRadiationInMat::GetIntProbSR(), and G4Integrator< T, F >::Laguerre().

00842 {
00843    const G4double tolerance = 1.0e-10 ;
00844    const G4int maxNumber = 12 ;
00845    G4int i, j, k ;
00846    G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ;
00847    G4double integral = 0.0 ;
00848 
00849    G4int fNumber = nLaguerre ;
00850    G4double* fAbscissa = new G4double[fNumber] ;
00851    G4double* fWeight   = new G4double[fNumber] ;
00852       
00853    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
00854    {
00855       if(i == 1)
00856       {
00857          nwt = (1.0 + alpha)*(3.0 + 0.92*alpha)
00858                 / (1.0 + 2.4*fNumber + 1.8*alpha) ;
00859       }
00860       else if(i == 2)
00861       {
00862          nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
00863       }
00864       else
00865       {
00866          cofi = i - 2 ;
00867          nwt += ((1.0+2.55*cofi)/(1.9*cofi)
00868               + 1.26*cofi*alpha/(1.0+3.5*cofi))
00869               * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
00870       }
00871       for(k=1;k<=maxNumber;k++)
00872       {
00873          temp1 = 1.0 ;
00874          temp2 = 0.0 ;
00875 
00876          for(j=1;j<=fNumber;j++)
00877          {
00878             temp3 = temp2 ;
00879             temp2 = temp1 ;
00880          temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ;
00881          }
00882          temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ;
00883          nwt1 = nwt ;
00884          nwt  = nwt1 - temp1/temp ;
00885 
00886          if(std::fabs(nwt - nwt1) <= tolerance) 
00887          {
00888             break ;
00889          }
00890       }
00891       if(k > maxNumber)
00892       {
00893          G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error",
00894                      FatalException, "Too many (>12) iterations.");
00895       }
00896          
00897       fAbscissa[i-1] =  nwt ;
00898       fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 
00899                 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
00900    }
00901 
00902    //
00903    // Integral evaluation
00904    //
00905 
00906    for(i=0;i<fNumber;i++)
00907    {
00908       integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ;
00909    }
00910    delete[] fAbscissa;
00911    delete[] fWeight;
00912    return integral ;
00913 }

template<class T, class F>
G4double G4Integrator< T, F >::Legendre ( G4double(*)(G4double f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 398 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4INCL::Math::pi.

00400 {
00401    G4double nwt, nwt1, temp1, temp2, temp3, temp ;
00402    G4double xDiff, xMean, dx, integral ;
00403 
00404    const G4double tolerance = 1.6e-10 ;
00405    G4int i, j,   k = nLegendre ;
00406    G4int fNumber = (nLegendre + 1)/2 ;
00407 
00408    if(2*fNumber != k)
00409    {
00410       G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall",
00411                   FatalException, "Invalid (odd) nLegendre in constructor.");
00412    }
00413 
00414    G4double* fAbscissa = new G4double[fNumber] ;
00415    G4double* fWeight   = new G4double[fNumber] ;
00416       
00417    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
00418    {
00419       nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ;  // Initial root approximation
00420 
00421       do     // loop of Newton's method  
00422       {                           
00423          temp1 = 1.0 ;
00424          temp2 = 0.0 ;
00425          for(j=1;j<=k;j++)
00426          {
00427             temp3 = temp2 ;
00428             temp2 = temp1 ;
00429             temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ;
00430          }
00431          temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ;
00432          nwt1 = nwt ;
00433          nwt  = nwt1 - temp1/temp ;       // Newton's method
00434       }
00435       while(std::fabs(nwt - nwt1) > tolerance) ;
00436          
00437       fAbscissa[fNumber-i] =  nwt ;
00438       fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ;
00439    }
00440 
00441    //
00442    // Now we ready to get integral 
00443    //
00444    
00445    xMean = 0.5*(a + b) ;
00446    xDiff = 0.5*(b - a) ;
00447    integral = 0.0 ;
00448    for(i=0;i<fNumber;i++)
00449    {
00450       dx = xDiff*fAbscissa[i] ;
00451       integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)    ) ;
00452    }
00453    delete[] fAbscissa;
00454    delete[] fWeight;
00455 
00456    return integral *= xDiff ;
00457 } 

template<class T, class F>
G4double G4Integrator< T, F >::Legendre ( T *  ptrT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 387 of file G4Integrator.icc.

References G4Integrator< T, F >::Legendre().

00389 {
00390   return Legendre(*ptrT,f,a,b,nLegendre) ;
00391 }

template<class T, class F>
G4double G4Integrator< T, F >::Legendre ( T &  typeT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 321 of file G4Integrator.icc.

References FatalException, G4Exception(), and G4INCL::Math::pi.

Referenced by G4Integrator< T, F >::Legendre().

00323 {
00324    G4double nwt, nwt1, temp1, temp2, temp3, temp ;
00325    G4double xDiff, xMean, dx, integral ;
00326 
00327    const G4double tolerance = 1.6e-10 ;
00328    G4int i, j,   k = nLegendre ;
00329    G4int fNumber = (nLegendre + 1)/2 ;
00330 
00331    if(2*fNumber != k)
00332    {
00333       G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall",
00334                   FatalException, "Invalid (odd) nLegendre in constructor.");
00335    }
00336 
00337    G4double* fAbscissa = new G4double[fNumber] ;
00338    G4double* fWeight   = new G4double[fNumber] ;
00339       
00340    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
00341    {
00342       nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ;  // Initial root approximation
00343 
00344       do     // loop of Newton's method  
00345       {                           
00346          temp1 = 1.0 ;
00347          temp2 = 0.0 ;
00348          for(j=1;j<=k;j++)
00349          {
00350             temp3 = temp2 ;
00351             temp2 = temp1 ;
00352             temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ;
00353          }
00354          temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ;
00355          nwt1 = nwt ;
00356          nwt  = nwt1 - temp1/temp ;       // Newton's method
00357       }
00358       while(std::fabs(nwt - nwt1) > tolerance) ;
00359          
00360       fAbscissa[fNumber-i] =  nwt ;
00361       fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ;
00362    }
00363 
00364    //
00365    // Now we ready to get integral 
00366    //
00367    
00368    xMean = 0.5*(a + b) ;
00369    xDiff = 0.5*(b - a) ;
00370    integral = 0.0 ;
00371    for(i=0;i<fNumber;i++)
00372    {
00373       dx = xDiff*fAbscissa[i] ;
00374       integral += fWeight[i]*( (typeT.*f)(xMean + dx) + 
00375                                (typeT.*f)(xMean - dx)    ) ;
00376    }
00377    delete[] fAbscissa;
00378    delete[] fWeight;
00379    return integral *= xDiff ;
00380 } 

template<class T, class F>
G4double G4Integrator< T, F >::Legendre10 ( G4double(*)(G4double f,
G4double  a,
G4double  b 
)

Definition at line 509 of file G4Integrator.icc.

00511 {
00512    G4int i ;
00513    G4double xDiff, xMean, dx, integral ;
00514    
00515    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
00516    
00517    static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
00518                                   0.679409568299024, 0.865063366688985,
00519                                   0.973906528517172                      } ;
00520    
00521    static G4double weight[] =   { 0.295524224714753, 0.269266719309996, 
00522                                   0.219086362515982, 0.149451349150581,
00523                                   0.066671344308688                      } ;
00524    xMean = 0.5*(a + b) ;
00525    xDiff = 0.5*(b - a) ;
00526    integral = 0.0 ;
00527    for(i=0;i<5;i++)
00528    {
00529      dx = xDiff*abscissa[i] ;
00530      integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ;
00531    }
00532    return integral *= xDiff ;
00533 }

template<class T, class F>
G4double G4Integrator< T, F >::Legendre10 ( T *  ptrT,
f,
G4double  a,
G4double  b 
)

Definition at line 499 of file G4Integrator.icc.

References G4Integrator< T, F >::Legendre10().

00500 {
00501   return Legendre10(*ptrT,f,a,b) ;
00502 } 

template<class T, class F>
G4double G4Integrator< T, F >::Legendre10 ( T &  typeT,
f,
G4double  a,
G4double  b 
)

Definition at line 469 of file G4Integrator.icc.

Referenced by G4VXTRenergyLoss::BuildAngleForEnergyBank(), G4NuclNuclDiffuseElastic::BuildAngleTable(), G4DiffuseElastic::BuildAngleTable(), G4VXTRenergyLoss::BuildEnergyTable(), G4InitXscPAI::IntegralCherenkov(), G4InitXscPAI::IntegralPAIdEdx(), G4InitXscPAI::IntegralPAIxSection(), G4InitXscPAI::IntegralPlasmon(), G4Integrator< T, F >::Legendre10(), G4NuclNuclDiffuseElastic::SampleThetaCMS(), G4DiffuseElastic::SampleThetaCMS(), G4NuclNuclDiffuseElastic::TestAngleTable(), G4DiffuseElastic::TestAngleTable(), and G4VXTRenergyLoss::XTRNSpectralDensity().

00470 {
00471    G4int i ;
00472    G4double xDiff, xMean, dx, integral ;
00473    
00474    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916
00475    
00476    static G4double abscissa[] = { 0.148874338981631, 0.433395394129247,
00477                                   0.679409568299024, 0.865063366688985,
00478                                   0.973906528517172                      } ;
00479    
00480    static G4double weight[] =   { 0.295524224714753, 0.269266719309996, 
00481                                   0.219086362515982, 0.149451349150581,
00482                                   0.066671344308688                      } ;
00483    xMean = 0.5*(a + b) ;
00484    xDiff = 0.5*(b - a) ;
00485    integral = 0.0 ;
00486    for(i=0;i<5;i++)
00487    {
00488      dx = xDiff*abscissa[i] ;
00489      integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
00490    }
00491    return integral *= xDiff ;
00492 }

template<class T, class F>
G4double G4Integrator< T, F >::Legendre96 ( G4double(*)(G4double f,
G4double  a,
G4double  b 
)

Definition at line 647 of file G4Integrator.icc.

00649 {
00650    G4int i ;
00651    G4double xDiff, xMean, dx, integral ;
00652    
00653    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
00654    
00655    static G4double 
00656    abscissa[] = { 
00657                   0.016276744849602969579, 0.048812985136049731112,
00658                   0.081297495464425558994, 0.113695850110665920911,
00659                   0.145973714654896941989, 0.178096882367618602759,  // 6
00660                            
00661                   0.210031310460567203603, 0.241743156163840012328,
00662                   0.273198812591049141487, 0.304364944354496353024,
00663                   0.335208522892625422616, 0.365696861472313635031,  // 12
00664                            
00665                   0.395797649828908603285, 0.425478988407300545365,
00666                   0.454709422167743008636, 0.483457973920596359768,
00667                   0.511694177154667673586, 0.539388108324357436227,  // 18
00668                            
00669                   0.566510418561397168404, 0.593032364777572080684,
00670                   0.618925840125468570386, 0.644163403784967106798,
00671                   0.668718310043916153953, 0.692564536642171561344,  // 24
00672                            
00673                   0.715676812348967626225, 0.738030643744400132851,
00674                   0.759602341176647498703, 0.780369043867433217604,
00675                   0.800308744139140817229, 0.819400310737931675539,  // 30
00676                            
00677                   0.837623511228187121494, 0.854959033434601455463,
00678                   0.871388505909296502874, 0.886894517402420416057,
00679                   0.901460635315852341319, 0.915071423120898074206,  // 36
00680                            
00681                   0.927712456722308690965, 0.939370339752755216932,
00682                   0.950032717784437635756, 0.959688291448742539300,
00683                   0.968326828463264212174, 0.975939174585136466453,  // 42
00684                            
00685                   0.982517263563014677447, 0.988054126329623799481,
00686                   0.992543900323762624572, 0.995981842987209290650,
00687                   0.998364375863181677724, 0.999689503883230766828   // 48
00688                                                                             } ;
00689    
00690    static G4double 
00691    weight[] = {  
00692                   0.032550614492363166242, 0.032516118713868835987,
00693                   0.032447163714064269364, 0.032343822568575928429,
00694                   0.032206204794030250669, 0.032034456231992663218,  // 6
00695                            
00696                   0.031828758894411006535, 0.031589330770727168558,
00697                   0.031316425596862355813, 0.031010332586313837423,
00698                   0.030671376123669149014, 0.030299915420827593794,  // 12
00699                            
00700                   0.029896344136328385984, 0.029461089958167905970,
00701                   0.028994614150555236543, 0.028497411065085385646,
00702                   0.027970007616848334440, 0.027412962726029242823,  // 18
00703                            
00704                   0.026826866725591762198, 0.026212340735672413913,
00705                   0.025570036005349361499, 0.024900633222483610288,
00706                   0.024204841792364691282, 0.023483399085926219842,  // 24
00707                            
00708                   0.022737069658329374001, 0.021966644438744349195,
00709                   0.021172939892191298988, 0.020356797154333324595,
00710                   0.019519081140145022410, 0.018660679627411467385,  // 30
00711                            
00712                   0.017782502316045260838, 0.016885479864245172450,
00713                   0.015970562902562291381, 0.015038721026994938006,
00714                   0.014090941772314860916, 0.013128229566961572637,  // 36
00715                            
00716                   0.012151604671088319635, 0.011162102099838498591,
00717                   0.010160770535008415758, 0.009148671230783386633,
00718                   0.008126876925698759217, 0.007096470791153865269,  // 42
00719                            
00720                   0.006058545504235961683, 0.005014202742927517693,
00721                   0.003964554338444686674, 0.002910731817934946408,
00722                   0.001853960788946921732, 0.000796792065552012429   // 48
00723                                                                             } ;
00724    xMean = 0.5*(a + b) ;
00725    xDiff = 0.5*(b - a) ;
00726    integral = 0.0 ;
00727    for(i=0;i<48;i++)
00728    {
00729       dx = xDiff*abscissa[i] ;
00730       integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ;
00731    }
00732    return integral *= xDiff ;
00733 }

template<class T, class F>
G4double G4Integrator< T, F >::Legendre96 ( T *  ptrT,
f,
G4double  a,
G4double  b 
)

Definition at line 637 of file G4Integrator.icc.

References G4Integrator< T, F >::Legendre96().

00638 {
00639   return Legendre96(*ptrT,f,a,b) ;
00640 } 

template<class T, class F>
G4double G4Integrator< T, F >::Legendre96 ( T &  typeT,
f,
G4double  a,
G4double  b 
)

Definition at line 545 of file G4Integrator.icc.

Referenced by G4VXTRenergyLoss::BuildGlobalAngleTable(), G4NuclNuclDiffuseElastic::GetCint(), G4NuclNuclDiffuseElastic::GetErfInt(), G4NuclNuclDiffuseElastic::GetSint(), G4NuclNuclDiffuseElastic::IntegralElasticProb(), G4DiffuseElastic::IntegralElasticProb(), G4Integrator< T, F >::Legendre96(), G4NuclNuclDiffuseElastic::SampleThetaCMS(), G4DiffuseElastic::SampleThetaCMS(), G4VXTRenergyLoss::SpectralXTRdEdx(), G4NuclNuclDiffuseElastic::TestAngleTable(), G4DiffuseElastic::TestAngleTable(), G4VXTRenergyLoss::XTRNAngleDensity(), and G4VXTRenergyLoss::XTRNSpectralDensity().

00546 {
00547    G4int i ;
00548    G4double xDiff, xMean, dx, integral ;
00549    
00550    // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919
00551    
00552    static G4double 
00553    abscissa[] = { 
00554                   0.016276744849602969579, 0.048812985136049731112,
00555                   0.081297495464425558994, 0.113695850110665920911,
00556                   0.145973714654896941989, 0.178096882367618602759,  // 6
00557                            
00558                   0.210031310460567203603, 0.241743156163840012328,
00559                   0.273198812591049141487, 0.304364944354496353024,
00560                   0.335208522892625422616, 0.365696861472313635031,  // 12
00561                            
00562                   0.395797649828908603285, 0.425478988407300545365,
00563                   0.454709422167743008636, 0.483457973920596359768,
00564                   0.511694177154667673586, 0.539388108324357436227,  // 18
00565                            
00566                   0.566510418561397168404, 0.593032364777572080684,
00567                   0.618925840125468570386, 0.644163403784967106798,
00568                   0.668718310043916153953, 0.692564536642171561344,  // 24
00569                            
00570                   0.715676812348967626225, 0.738030643744400132851,
00571                   0.759602341176647498703, 0.780369043867433217604,
00572                   0.800308744139140817229, 0.819400310737931675539,  // 30
00573                            
00574                   0.837623511228187121494, 0.854959033434601455463,
00575                   0.871388505909296502874, 0.886894517402420416057,
00576                   0.901460635315852341319, 0.915071423120898074206,  // 36
00577                            
00578                   0.927712456722308690965, 0.939370339752755216932,
00579                   0.950032717784437635756, 0.959688291448742539300,
00580                   0.968326828463264212174, 0.975939174585136466453,  // 42
00581                            
00582                   0.982517263563014677447, 0.988054126329623799481,
00583                   0.992543900323762624572, 0.995981842987209290650,
00584                   0.998364375863181677724, 0.999689503883230766828   // 48
00585                                                                             } ;
00586    
00587    static G4double 
00588    weight[] = {  
00589                   0.032550614492363166242, 0.032516118713868835987,
00590                   0.032447163714064269364, 0.032343822568575928429,
00591                   0.032206204794030250669, 0.032034456231992663218,  // 6
00592                            
00593                   0.031828758894411006535, 0.031589330770727168558,
00594                   0.031316425596862355813, 0.031010332586313837423,
00595                   0.030671376123669149014, 0.030299915420827593794,  // 12
00596                            
00597                   0.029896344136328385984, 0.029461089958167905970,
00598                   0.028994614150555236543, 0.028497411065085385646,
00599                   0.027970007616848334440, 0.027412962726029242823,  // 18
00600                            
00601                   0.026826866725591762198, 0.026212340735672413913,
00602                   0.025570036005349361499, 0.024900633222483610288,
00603                   0.024204841792364691282, 0.023483399085926219842,  // 24
00604                            
00605                   0.022737069658329374001, 0.021966644438744349195,
00606                   0.021172939892191298988, 0.020356797154333324595,
00607                   0.019519081140145022410, 0.018660679627411467385,  // 30
00608                            
00609                   0.017782502316045260838, 0.016885479864245172450,
00610                   0.015970562902562291381, 0.015038721026994938006,
00611                   0.014090941772314860916, 0.013128229566961572637,  // 36
00612                            
00613                   0.012151604671088319635, 0.011162102099838498591,
00614                   0.010160770535008415758, 0.009148671230783386633,
00615                   0.008126876925698759217, 0.007096470791153865269,  // 42
00616                            
00617                   0.006058545504235961683, 0.005014202742927517693,
00618                   0.003964554338444686674, 0.002910731817934946408,
00619                   0.001853960788946921732, 0.000796792065552012429   // 48
00620                                                                             } ;
00621    xMean = 0.5*(a + b) ;
00622    xDiff = 0.5*(b - a) ;
00623    integral = 0.0 ;
00624    for(i=0;i<48;i++)
00625    {
00626       dx = xDiff*abscissa[i] ;
00627       integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ;
00628    }
00629    return integral *= xDiff ;
00630 }

template<class T, class F>
G4double G4Integrator< T, F >::Simpson ( G4double(*)(G4double f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 105 of file G4Integrator.icc.

00109 {
00110    G4int    i ;
00111    G4double step = (xFinal - xInitial)/iterationNumber ;
00112    G4double x = xInitial ;
00113    G4double xPlus = xInitial + 0.5*step ;
00114    G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ;
00115    G4double sum = (*f)(xPlus) ;
00116 
00117    for(i=1;i<iterationNumber;i++)
00118    {
00119       x     += step ;
00120       xPlus += step ;
00121       mean  += (*f)(x) ;
00122       sum   += (*f)(xPlus) ;
00123    }
00124    mean += 2.0*sum ;
00125 
00126    return mean*step/3.0 ;   
00127 }

template<class T, class F>
G4double G4Integrator< T, F >::Simpson ( T *  ptrT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 73 of file G4Integrator.icc.

00078 {
00079    G4int    i ;
00080    G4double step = (xFinal - xInitial)/iterationNumber ;
00081    G4double x = xInitial ;
00082    G4double xPlus = xInitial + 0.5*step ;
00083    G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ;
00084    G4double sum = (ptrT->*f)(xPlus) ;
00085 
00086    for(i=1;i<iterationNumber;i++)
00087    {
00088       x     += step ;
00089       xPlus += step ;
00090       mean  += (ptrT->*f)(x) ;
00091       sum   += (ptrT->*f)(xPlus) ;
00092    }
00093    mean += 2.0*sum ;
00094 
00095    return mean*step/3.0 ;   
00096 }

template<class T, class F>
G4double G4Integrator< T, F >::Simpson ( T &  typeT,
f,
G4double  a,
G4double  b,
G4int  n 
)

Definition at line 42 of file G4Integrator.icc.

Referenced by G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(), G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(), G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(), and G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond().

00047 {
00048    G4int    i ;
00049    G4double step = (xFinal - xInitial)/iterationNumber ;
00050    G4double x = xInitial ;
00051    G4double xPlus = xInitial + 0.5*step ;
00052    G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ;
00053    G4double sum = (typeT.*f)(xPlus) ;
00054 
00055    for(i=1;i<iterationNumber;i++)
00056    {
00057       x     += step ;
00058       xPlus += step ;
00059       mean  += (typeT.*f)(x) ;
00060       sum   += (typeT.*f)(xPlus) ;
00061    }
00062    mean += 2.0*sum ;
00063 
00064    return mean*step/3.0 ;   
00065 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:16 2013 for Geant4 by  doxygen 1.4.7