G4GaussHermiteQ Class Reference

#include <G4GaussHermiteQ.hh>

Inheritance diagram for G4GaussHermiteQ:

G4VGaussianQuadrature

Public Member Functions

 G4GaussHermiteQ (function pFunction, G4int nHermite)
G4double Integral () const

Detailed Description

Definition at line 60 of file G4GaussHermiteQ.hh.


Constructor & Destructor Documentation

G4GaussHermiteQ::G4GaussHermiteQ ( function  pFunction,
G4int  nHermite 
)

Definition at line 37 of file G4GaussHermiteQ.cc.

References G4VGaussianQuadrature::fAbscissa, FatalException, G4VGaussianQuadrature::fNumber, G4VGaussianQuadrature::fWeight, G4Exception(), and G4INCL::Math::pi.

00039    : G4VGaussianQuadrature(pFunction)
00040 {
00041    const G4double tolerance = 1.0e-12 ;
00042    const G4int maxNumber = 12 ;
00043    
00044    G4int i=1, j=1, k=1 ;
00045    G4double newton0=0.;
00046    G4double newton1=0.0, temp1=0.0, temp2=0.0, temp3=0.0, temp=0.0 ;
00047    G4double piInMinusQ = std::pow(pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ??
00048 
00049    fNumber = (nHermite +1)/2 ;
00050    fAbscissa = new G4double[fNumber] ;
00051    fWeight   = new G4double[fNumber] ;
00052 
00053    for(i=1;i<=fNumber;i++)
00054    {
00055       if(i == 1)
00056       {
00057          newton0 = std::sqrt((G4double)(2*nHermite + 1)) - 
00058                    1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ;
00059       }
00060       else if(i == 2)
00061       {
00062          newton0 -= 1.14001*std::pow((G4double)nHermite,0.425999)/newton0 ;
00063       }
00064       else if(i == 3)
00065       {
00066          newton0 = 1.86002*newton0 - 0.86002*fAbscissa[0] ;
00067       }
00068       else if(i == 4)
00069       {
00070          newton0 = 1.91001*newton0 - 0.91001*fAbscissa[1] ;
00071       }
00072       else 
00073       {
00074          newton0 = 2.0*newton0 - fAbscissa[i - 3] ;
00075       }
00076       for(k=1;k<=maxNumber;k++)
00077       {
00078          temp1 = piInMinusQ ;
00079          temp2 = 0.0 ;
00080          for(j=1;j<=nHermite;j++)
00081          {
00082             temp3 = temp2 ;
00083             temp2 = temp1 ;
00084             temp1 = newton0*std::sqrt(2.0/j)*temp2
00085                   - std::sqrt(((G4double)(j - 1))/j)*temp3 ;
00086          }
00087          temp = std::sqrt((G4double)2*nHermite)*temp2 ;
00088          newton1 = newton0 ;
00089          newton0 = newton1 - temp1/temp ;
00090          if(std::fabs(newton0 - newton1) <= tolerance) 
00091          {
00092             break ;
00093          }
00094       }
00095       if(k > maxNumber)
00096       {
00097          G4Exception("G4GaussHermiteQ::G4GaussHermiteQ()",
00098                      "OutOfRange", FatalException,
00099                      "Too many iterations in Gauss-Hermite constructor.") ;
00100       }
00101       fAbscissa[i-1] =  newton0 ;
00102       fWeight[i-1] = 2.0/(temp*temp) ;
00103    }
00104 }


Member Function Documentation

G4double G4GaussHermiteQ::Integral (  )  const

Definition at line 112 of file G4GaussHermiteQ.cc.

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

00113 {
00114    G4double integral = 0.0 ;
00115    for(G4int i=0;i<fNumber;i++)
00116    {
00117       integral += fWeight[i]*(fFunction(fAbscissa[i])
00118                 + fFunction(-fAbscissa[i])) ;
00119    }
00120    return integral ;
00121 }


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