00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4GaussLaguerreQ.hh"
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 G4GaussLaguerreQ::G4GaussLaguerreQ( function pFunction,
00043 G4double alpha,
00044 G4int nLaguerre )
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++)
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 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 G4double
00116 G4GaussLaguerreQ::Integral() const
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 }