#include <G4GaussJacobiQ.hh>
Inheritance diagram for G4GaussJacobiQ:
Public Member Functions | |
G4GaussJacobiQ (function pFunction, G4double alpha, G4double beta, G4int nJacobi) | |
G4double | Integral () const |
Definition at line 62 of file G4GaussJacobiQ.hh.
Definition at line 37 of file G4GaussJacobiQ.cc.
References G4VGaussianQuadrature::fAbscissa, FatalException, G4VGaussianQuadrature::fNumber, G4VGaussianQuadrature::fWeight, G4Exception(), and G4VGaussianQuadrature::GammaLogarithm().
00041 : G4VGaussianQuadrature(pFunction) 00042 00043 { 00044 const G4double tolerance = 1.0e-12 ; 00045 const G4double maxNumber = 12 ; 00046 G4int i=1, k=1 ; 00047 G4double root=0.; 00048 G4double alphaBeta=0.0, alphaReduced=0.0, betaReduced=0.0, 00049 root1=0.0, root2=0.0, root3=0.0 ; 00050 G4double a=0.0, b=0.0, c=0.0, 00051 newton1=0.0, newton2=0.0, newton3=0.0, newton0=0.0, 00052 temp=0.0, rootTemp=0.0 ; 00053 00054 fNumber = nJacobi ; 00055 fAbscissa = new G4double[fNumber] ; 00056 fWeight = new G4double[fNumber] ; 00057 00058 for (i=1;i<=nJacobi;i++) 00059 { 00060 if (i == 1) 00061 { 00062 alphaReduced = alpha/nJacobi ; 00063 betaReduced = beta/nJacobi ; 00064 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+ 00065 0.767999*alphaReduced/nJacobi) ; 00066 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced 00067 + 0.451998*alphaReduced*alphaReduced 00068 + 0.83001*alphaReduced*betaReduced ; 00069 root = 1.0-root1/root2 ; 00070 } 00071 else if (i == 2) 00072 { 00073 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ; 00074 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ; 00075 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ; 00076 root -= (1.0-root)*root1*root2*root3 ; 00077 } 00078 else if (i == 3) 00079 { 00080 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ; 00081 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ; 00082 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ; 00083 root -= (fAbscissa[0]-root)*root1*root2*root3 ; 00084 } 00085 else if (i == nJacobi-1) 00086 { 00087 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ; 00088 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ; 00089 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ; 00090 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ; 00091 } 00092 else if (i == nJacobi) 00093 { 00094 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ; 00095 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ; 00096 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ; 00097 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ; 00098 } 00099 else 00100 { 00101 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ; 00102 } 00103 alphaBeta = alpha + beta ; 00104 for (k=1;k<=maxNumber;k++) 00105 { 00106 temp = 2.0 + alphaBeta ; 00107 newton1 = (alpha-beta+temp*root)/2.0 ; 00108 newton2 = 1.0 ; 00109 for (G4int j=2;j<=nJacobi;j++) 00110 { 00111 newton3 = newton2 ; 00112 newton2 = newton1 ; 00113 temp = 2*j+alphaBeta ; 00114 a = 2*j*(j+alphaBeta)*(temp-2.0) ; 00115 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ; 00116 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ; 00117 newton1 = (b*newton2-c*newton3)/a ; 00118 } 00119 newton0 = (nJacobi*(alpha - beta - temp*root)*newton1 + 00120 2.0*(nJacobi + alpha)*(nJacobi + beta)*newton2)/ 00121 (temp*(1.0 - root*root)) ; 00122 rootTemp = root ; 00123 root = rootTemp - newton1/newton0 ; 00124 if (std::fabs(root-rootTemp) <= tolerance) 00125 { 00126 break ; 00127 } 00128 } 00129 if (k > maxNumber) 00130 { 00131 G4Exception("G4GaussJacobiQ::G4GaussJacobiQ()", "OutOfRange", 00132 FatalException, "Too many iterations in constructor.") ; 00133 } 00134 fAbscissa[i-1] = root ; 00135 fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 00136 GammaLogarithm((G4double)(beta+nJacobi)) - 00137 GammaLogarithm((G4double)(nJacobi+1.0)) - 00138 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) 00139 *temp*std::pow(2.0,alphaBeta)/(newton0*newton2) ; 00140 } 00141 }
G4double G4GaussJacobiQ::Integral | ( | ) | const |
Definition at line 152 of file G4GaussJacobiQ.cc.
References G4VGaussianQuadrature::fAbscissa, G4VGaussianQuadrature::fFunction, G4VGaussianQuadrature::fNumber, and G4VGaussianQuadrature::fWeight.
00153 { 00154 G4double integral = 0.0 ; 00155 for(G4int i=0;i<fNumber;i++) 00156 { 00157 integral += fWeight[i]*fFunction(fAbscissa[i]) ; 00158 } 00159 return integral ; 00160 }