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 "G4GaussJacobiQ.hh"
00030
00031
00032
00033
00034
00035
00036
00037 G4GaussJacobiQ::G4GaussJacobiQ( function pFunction,
00038 G4double alpha,
00039 G4double beta,
00040 G4int nJacobi )
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 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151 G4double
00152 G4GaussJacobiQ::Integral() const
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 }
00161