#include <G4GaussLaguerreQ.hh>
Inheritance diagram for G4GaussLaguerreQ:
Public Member Functions | |
G4GaussLaguerreQ (function pFunction, G4double alpha, G4int nLaguerre) | |
G4double | Integral () const |
Definition at line 68 of file G4GaussLaguerreQ.hh.
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 }
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 }