G4GaussLaguerreQ Class Reference

#include <G4GaussLaguerreQ.hh>

Inheritance diagram for G4GaussLaguerreQ:

G4VGaussianQuadrature

Public Member Functions

 G4GaussLaguerreQ (function pFunction, G4double alpha, G4int nLaguerre)
G4double Integral () const

Detailed Description

Definition at line 68 of file G4GaussLaguerreQ.hh.


Constructor & Destructor Documentation

G4GaussLaguerreQ::G4GaussLaguerreQ ( function  pFunction,
G4double  alpha,
G4int  nLaguerre 
)

Definition at line 42 of file G4GaussLaguerreQ.cc.

References G4VGaussianQuadrature::fAbscissa, FatalException, G4VGaussianQuadrature::fNumber, G4VGaussianQuadrature::fWeight, G4Exception(), and G4VGaussianQuadrature::GammaLogarithm().

00045    : G4VGaussianQuadrature(pFunction)
00046 {
00047    const G4double tolerance = 1.0e-10 ;
00048    const G4int maxNumber = 12 ;
00049    G4int i=1, k=1 ;
00050    G4double newton0=0.0, newton1=0.0,
00051             temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0, cofi=0.0 ;
00052 
00053    fNumber = nLaguerre ;
00054    fAbscissa = new G4double[fNumber] ;
00055    fWeight   = new G4double[fNumber] ;
00056       
00057    for(i=1;i<=fNumber;i++)      // Loop over the desired roots
00058    {
00059       if(i == 1)
00060       {
00061          newton0 = (1.0 + alpha)*(3.0 + 0.92*alpha)
00062                  / (1.0 + 2.4*fNumber + 1.8*alpha) ;
00063       }
00064       else if(i == 2)
00065       {
00066          newton0 += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ;
00067       }
00068       else
00069       {
00070          cofi = i - 2 ;
00071          newton0 += ((1.0+2.55*cofi)/(1.9*cofi)
00072                     + 1.26*cofi*alpha/(1.0+3.5*cofi))
00073                     * (newton0 - fAbscissa[i-3])/(1.0 + 0.3*alpha) ;
00074       }
00075       for(k=1;k<=maxNumber;k++)
00076       {
00077          temp1 = 1.0 ;
00078          temp2 = 0.0 ;
00079          for(G4int j=1;j<=fNumber;j++)
00080          {
00081             temp3 = temp2 ;
00082             temp2 = temp1 ;
00083             temp1 = ((2*j - 1 + alpha - newton0)*temp2
00084                      - (j - 1 + alpha)*temp3)/j ;
00085          }
00086          temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/newton0 ;
00087          newton1 = newton0 ;
00088          newton0 = newton1 - temp1/temp ;
00089          if(std::fabs(newton0 - newton1) <= tolerance) 
00090          {
00091             break ;
00092          }
00093       }
00094       if(k > maxNumber)
00095       {
00096          G4Exception("G4GaussLaguerreQ::G4GaussLaguerreQ()",
00097                      "OutOfRange", FatalException,
00098                      "Too many iterations in Gauss-Laguerre constructor") ;
00099       }
00100          
00101       fAbscissa[i-1] = newton0 ;
00102       fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber)
00103                    - GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ;
00104    }
00105 }


Member Function Documentation

G4double G4GaussLaguerreQ::Integral (  )  const

Definition at line 116 of file G4GaussLaguerreQ.cc.

References G4VGaussianQuadrature::fAbscissa, G4VGaussianQuadrature::fFunction, G4VGaussianQuadrature::fNumber, and G4VGaussianQuadrature::fWeight.

00117 {
00118    G4double integral = 0.0 ;
00119    for(G4int i=0;i<fNumber;i++)
00120    {
00121       integral += fWeight[i]*fFunction(fAbscissa[i]) ;
00122    }
00123    return integral ;
00124 }


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