#include <G4Integrator.hh>
Definition at line 49 of file G4Integrator.hh.
G4Integrator< T, F >::G4Integrator | ( | ) | [inline] |
G4Integrator< T, F >::~G4Integrator | ( | ) | [inline] |
void G4Integrator< T, F >::AdaptGauss | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b, | |||
G4double | e, | |||
G4double & | sum, | |||
G4int & | n | |||
) | [protected] |
Definition at line 227 of file G4Integrator.icc.
References G4Integrator< T, F >::AdaptGauss(), G4cout, G4endl, and G4Integrator< T, F >::Gauss().
00231 { 00232 if(depth > 100) 00233 { 00234 G4cout<<"G4SimpleIntegration::AdaptGauss: WARNING !!!"<<G4endl ; 00235 G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps " 00236 <<G4endl ; 00237 00238 return ; 00239 } 00240 G4double xMean = (xInitial + xFinal)/2.0 ; 00241 G4double leftHalf = Gauss(f,xInitial,xMean) ; 00242 G4double rightHalf = Gauss(f,xMean,xFinal) ; 00243 G4double full = Gauss(f,xInitial,xFinal) ; 00244 if(std::fabs(leftHalf+rightHalf-full) < fTolerance) 00245 { 00246 sum += full ; 00247 } 00248 else 00249 { 00250 depth++ ; 00251 AdaptGauss(f,xInitial,xMean,fTolerance,sum,depth) ; 00252 AdaptGauss(f,xMean,xFinal,fTolerance,sum,depth) ; 00253 } 00254 }
void G4Integrator< T, F >::AdaptGauss | ( | T * | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4double | e, | |||
G4double & | sum, | |||
G4int & | n | |||
) | [protected] |
Definition at line 215 of file G4Integrator.icc.
References G4Integrator< T, F >::AdaptGauss().
00219 { 00220 AdaptGauss(*ptrT,f,xInitial,xFinal,fTolerance,sum,depth) ; 00221 }
void G4Integrator< T, F >::AdaptGauss | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4double | e, | |||
G4double & | sum, | |||
G4int & | n | |||
) | [protected] |
Definition at line 185 of file G4Integrator.icc.
References G4cout, G4endl, and G4Integrator< T, F >::Gauss().
Referenced by G4Integrator< T, F >::AdaptGauss(), and G4Integrator< T, F >::AdaptiveGauss().
00189 { 00190 if(depth > 100) 00191 { 00192 G4cout<<"G4Integrator<T,F>::AdaptGauss: WARNING !!!"<<G4endl ; 00193 G4cout<<"Function varies too rapidly to get stated accuracy in 100 steps " 00194 <<G4endl ; 00195 00196 return ; 00197 } 00198 G4double xMean = (xInitial + xFinal)/2.0 ; 00199 G4double leftHalf = Gauss(typeT,f,xInitial,xMean) ; 00200 G4double rightHalf = Gauss(typeT,f,xMean,xFinal) ; 00201 G4double full = Gauss(typeT,f,xInitial,xFinal) ; 00202 if(std::fabs(leftHalf+rightHalf-full) < fTolerance) 00203 { 00204 sum += full ; 00205 } 00206 else 00207 { 00208 depth++ ; 00209 AdaptGauss(typeT,f,xInitial,xMean,fTolerance,sum,depth) ; 00210 AdaptGauss(typeT,f,xMean,xFinal,fTolerance,sum,depth) ; 00211 } 00212 }
G4double G4Integrator< T, F >::AdaptiveGauss | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b, | |||
G4double | e | |||
) |
Definition at line 289 of file G4Integrator.icc.
References G4Integrator< T, F >::AdaptGauss().
00291 { 00292 G4int depth = 0 ; 00293 G4double sum = 0.0 ; 00294 AdaptGauss(f,xInitial,xFinal,e,sum,depth) ; 00295 return sum ; 00296 }
G4double G4Integrator< T, F >::AdaptiveGauss | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4double | e | |||
) |
Definition at line 277 of file G4Integrator.icc.
References G4Integrator< T, F >::AdaptiveGauss().
00279 { 00280 return AdaptiveGauss(*ptrT,f,xInitial,xFinal,e) ; 00281 }
G4double G4Integrator< T, F >::AdaptiveGauss | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4double | e | |||
) |
Definition at line 262 of file G4Integrator.icc.
References G4Integrator< T, F >::AdaptGauss().
Referenced by G4Integrator< T, F >::AdaptiveGauss(), G4NuclNuclDiffuseElastic::TestAngleTable(), and G4DiffuseElastic::TestAngleTable().
00264 { 00265 G4int depth = 0 ; 00266 G4double sum = 0.0 ; 00267 AdaptGauss(typeT,f,xInitial,xFinal,e,sum,depth) ; 00268 return sum ; 00269 }
G4double G4Integrator< T, F >::Chebyshev | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 794 of file G4Integrator.icc.
References G4INCL::Math::pi.
00796 { 00797 G4int i ; 00798 G4double xDiff, xMean, dx, integral = 0.0 ; 00799 00800 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ?? 00801 G4double cof = CLHEP::pi/fNumber ; 00802 G4double* fAbscissa = new G4double[fNumber] ; 00803 G4double* fWeight = new G4double[fNumber] ; 00804 for(i=0;i<fNumber;i++) 00805 { 00806 fAbscissa[i] = std::cos(cof*(i + 0.5)) ; 00807 fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ; 00808 } 00809 00810 // 00811 // Now we ready to estimate the integral 00812 // 00813 00814 xMean = 0.5*(a + b) ; 00815 xDiff = 0.5*(b - a) ; 00816 for(i=0;i<fNumber;i++) 00817 { 00818 dx = xDiff*fAbscissa[i] ; 00819 integral += fWeight[i]*(*f)(xMean + dx) ; 00820 } 00821 delete[] fAbscissa; 00822 delete[] fWeight; 00823 return integral *= xDiff ; 00824 }
G4double G4Integrator< T, F >::Chebyshev | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
G4double G4Integrator< T, F >::Chebyshev | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 746 of file G4Integrator.icc.
References G4INCL::Math::pi.
Referenced by G4Integrator< T, F >::Chebyshev().
00748 { 00749 G4int i ; 00750 G4double xDiff, xMean, dx, integral = 0.0 ; 00751 00752 G4int fNumber = nChebyshev ; // Try to reduce fNumber twice ?? 00753 G4double cof = CLHEP::pi/fNumber ; 00754 G4double* fAbscissa = new G4double[fNumber] ; 00755 G4double* fWeight = new G4double[fNumber] ; 00756 for(i=0;i<fNumber;i++) 00757 { 00758 fAbscissa[i] = std::cos(cof*(i + 0.5)) ; 00759 fWeight[i] = cof*std::sqrt(1 - fAbscissa[i]*fAbscissa[i]) ; 00760 } 00761 00762 // 00763 // Now we ready to estimate the integral 00764 // 00765 00766 xMean = 0.5*(a + b) ; 00767 xDiff = 0.5*(b - a) ; 00768 for(i=0;i<fNumber;i++) 00769 { 00770 dx = xDiff*fAbscissa[i] ; 00771 integral += fWeight[i]*(typeT.*f)(xMean + dx) ; 00772 } 00773 delete[] fAbscissa; 00774 delete[] fWeight; 00775 return integral *= xDiff ; 00776 }
G4double G4Integrator< T, F >::GammaLogarithm | ( | G4double | xx | ) | [protected] |
Definition at line 1016 of file G4Integrator.icc.
Referenced by G4Integrator< T, F >::Jacobi(), and G4Integrator< T, F >::Laguerre().
01017 { 01018 static G4double cof[6] = { 76.18009172947146, -86.50532032941677, 01019 24.01409824083091, -1.231739572450155, 01020 0.1208650973866179e-2, -0.5395239384953e-5 } ; 01021 register G4int j; 01022 G4double x = xx - 1.0 ; 01023 G4double tmp = x + 5.5 ; 01024 tmp -= (x + 0.5) * std::log(tmp) ; 01025 G4double ser = 1.000000000190015 ; 01026 01027 for ( j = 0; j <= 5; j++ ) 01028 { 01029 x += 1.0 ; 01030 ser += cof[j]/x ; 01031 } 01032 return -tmp + std::log(2.5066282746310005*ser) ; 01033 }
G4double G4Integrator< T, F >::Gauss | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b | |||
) | [protected] |
Definition at line 167 of file G4Integrator.icc.
00169 { 00170 static G4double root = 1.0/std::sqrt(3.0) ; 00171 00172 G4double xMean = (xInitial + xFinal)/2.0 ; 00173 G4double Step = (xFinal - xInitial)/2.0 ; 00174 G4double delta = Step*root ; 00175 G4double sum = ( (*f)(xMean + delta) + (*f)(xMean - delta) ) ; 00176 00177 return sum*Step ; 00178 }
G4double G4Integrator< T, F >::Gauss | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) | [protected] |
Definition at line 157 of file G4Integrator.icc.
References G4Integrator< T, F >::Gauss().
00158 { 00159 return Gauss(*ptrT,f,a,b) ; 00160 }
G4double G4Integrator< T, F >::Gauss | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) | [protected] |
Definition at line 138 of file G4Integrator.icc.
Referenced by G4Integrator< T, F >::AdaptGauss(), and G4Integrator< T, F >::Gauss().
00140 { 00141 static G4double root = 1.0/std::sqrt(3.0) ; 00142 00143 G4double xMean = (xInitial + xFinal)/2.0 ; 00144 G4double Step = (xFinal - xInitial)/2.0 ; 00145 G4double delta = Step*root ; 00146 G4double sum = ((typeT.*f)(xMean + delta) + 00147 (typeT.*f)(xMean - delta)) ; 00148 00149 return sum*Step ; 00150 }
G4double G4Integrator< T, F >::Hermite | ( | G4double(*)(G4double) | f, | |
G4int | n | |||
) |
Definition at line 1145 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4INCL::Math::pi.
01146 { 01147 const G4double tolerance = 1.0e-12 ; 01148 const G4int maxNumber = 12 ; 01149 01150 G4int i, j, k ; 01151 G4double integral = 0.0 ; 01152 G4double nwt=0., nwt1, temp1, temp2, temp3, temp ; 01153 01154 G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ?? 01155 01156 G4int fNumber = (nHermite +1)/2 ; 01157 G4double* fAbscissa = new G4double[fNumber] ; 01158 G4double* fWeight = new G4double[fNumber] ; 01159 01160 for(i=1;i<=fNumber;i++) 01161 { 01162 if(i == 1) 01163 { 01164 nwt = std::sqrt((G4double)(2*nHermite + 1)) - 01165 1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ; 01166 } 01167 else if(i == 2) 01168 { 01169 nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ; 01170 } 01171 else if(i == 3) 01172 { 01173 nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ; 01174 } 01175 else if(i == 4) 01176 { 01177 nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ; 01178 } 01179 else 01180 { 01181 nwt = 2.0*nwt - fAbscissa[i - 3] ; 01182 } 01183 for(k=1;k<=maxNumber;k++) 01184 { 01185 temp1 = piInMinusQ ; 01186 temp2 = 0.0 ; 01187 01188 for(j=1;j<=nHermite;j++) 01189 { 01190 temp3 = temp2 ; 01191 temp2 = temp1 ; 01192 temp1 = nwt*std::sqrt(2.0/j)*temp2 - 01193 std::sqrt(((G4double)(j - 1))/j)*temp3 ; 01194 } 01195 temp = std::sqrt((G4double)2*nHermite)*temp2 ; 01196 nwt1 = nwt ; 01197 nwt = nwt1 - temp1/temp ; 01198 01199 if(std::fabs(nwt - nwt1) <= tolerance) 01200 { 01201 break ; 01202 } 01203 } 01204 if(k > maxNumber) 01205 { 01206 G4Exception("G4Integrator<T,F>::Hermite(...)", "Error", 01207 FatalException, "Too many (>12) iterations."); 01208 } 01209 fAbscissa[i-1] = nwt ; 01210 fWeight[i-1] = 2.0/(temp*temp) ; 01211 } 01212 01213 // 01214 // Integral calculation 01215 // 01216 01217 for(i=0;i<fNumber;i++) 01218 { 01219 integral += fWeight[i]*( (*f)(fAbscissa[i]) + (*f)(-fAbscissa[i]) ) ; 01220 } 01221 delete[] fAbscissa; 01222 delete[] fWeight; 01223 return integral ; 01224 }
G4double G4Integrator< T, F >::Hermite | ( | T * | ptrT, | |
F | f, | |||
G4int | n | |||
) |
G4double G4Integrator< T, F >::Hermite | ( | T & | typeT, | |
F | f, | |||
G4int | n | |||
) |
Definition at line 1047 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4INCL::Math::pi.
Referenced by G4Integrator< T, F >::Hermite().
01048 { 01049 const G4double tolerance = 1.0e-12 ; 01050 const G4int maxNumber = 12 ; 01051 01052 G4int i, j, k ; 01053 G4double integral = 0.0 ; 01054 G4double nwt=0., nwt1, temp1, temp2, temp3, temp ; 01055 01056 G4double piInMinusQ = std::pow(CLHEP::pi,-0.25) ; // 1.0/std::sqrt(std::sqrt(pi)) ?? 01057 01058 G4int fNumber = (nHermite +1)/2 ; 01059 G4double* fAbscissa = new G4double[fNumber] ; 01060 G4double* fWeight = new G4double[fNumber] ; 01061 01062 for(i=1;i<=fNumber;i++) 01063 { 01064 if(i == 1) 01065 { 01066 nwt = std::sqrt((G4double)(2*nHermite + 1)) - 01067 1.85575001*std::pow((G4double)(2*nHermite + 1),-0.16666999) ; 01068 } 01069 else if(i == 2) 01070 { 01071 nwt -= 1.14001*std::pow((G4double)nHermite,0.425999)/nwt ; 01072 } 01073 else if(i == 3) 01074 { 01075 nwt = 1.86002*nwt - 0.86002*fAbscissa[0] ; 01076 } 01077 else if(i == 4) 01078 { 01079 nwt = 1.91001*nwt - 0.91001*fAbscissa[1] ; 01080 } 01081 else 01082 { 01083 nwt = 2.0*nwt - fAbscissa[i - 3] ; 01084 } 01085 for(k=1;k<=maxNumber;k++) 01086 { 01087 temp1 = piInMinusQ ; 01088 temp2 = 0.0 ; 01089 01090 for(j=1;j<=nHermite;j++) 01091 { 01092 temp3 = temp2 ; 01093 temp2 = temp1 ; 01094 temp1 = nwt*std::sqrt(2.0/j)*temp2 - 01095 std::sqrt(((G4double)(j - 1))/j)*temp3 ; 01096 } 01097 temp = std::sqrt((G4double)2*nHermite)*temp2 ; 01098 nwt1 = nwt ; 01099 nwt = nwt1 - temp1/temp ; 01100 01101 if(std::fabs(nwt - nwt1) <= tolerance) 01102 { 01103 break ; 01104 } 01105 } 01106 if(k > maxNumber) 01107 { 01108 G4Exception("G4Integrator<T,F>::Hermite(T,F, ...)", "Error", 01109 FatalException, "Too many (>12) iterations."); 01110 } 01111 fAbscissa[i-1] = nwt ; 01112 fWeight[i-1] = 2.0/(temp*temp) ; 01113 } 01114 01115 // 01116 // Integral calculation 01117 // 01118 01119 for(i=0;i<fNumber;i++) 01120 { 01121 integral += fWeight[i]*( (typeT.*f)(fAbscissa[i]) + 01122 (typeT.*f)(-fAbscissa[i]) ) ; 01123 } 01124 delete[] fAbscissa; 01125 delete[] fWeight; 01126 return integral ; 01127 }
G4double G4Integrator< T, F >::Jacobi | ( | G4double(*)(G4double) | f, | |
G4double | alpha, | |||
G4double | beta, | |||
G4int | n | |||
) |
Definition at line 1365 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().
01367 { 01368 const G4double tolerance = 1.0e-12 ; 01369 const G4double maxNumber = 12 ; 01370 G4int i, k, j ; 01371 G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ; 01372 G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ; 01373 01374 G4int fNumber = nJacobi ; 01375 G4double* fAbscissa = new G4double[fNumber] ; 01376 G4double* fWeight = new G4double[fNumber] ; 01377 01378 for (i=1;i<=nJacobi;i++) 01379 { 01380 if (i == 1) 01381 { 01382 alphaReduced = alpha/nJacobi ; 01383 betaReduced = beta/nJacobi ; 01384 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+ 01385 0.767999*alphaReduced/nJacobi) ; 01386 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced + 01387 0.451998*alphaReduced*alphaReduced + 01388 0.83001*alphaReduced*betaReduced ; 01389 root = 1.0-root1/root2 ; 01390 } 01391 else if (i == 2) 01392 { 01393 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ; 01394 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ; 01395 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ; 01396 root -= (1.0-root)*root1*root2*root3 ; 01397 } 01398 else if (i == 3) 01399 { 01400 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ; 01401 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ; 01402 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ; 01403 root -= (fAbscissa[0]-root)*root1*root2*root3 ; 01404 } 01405 else if (i == nJacobi-1) 01406 { 01407 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ; 01408 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ; 01409 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ; 01410 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ; 01411 } 01412 else if (i == nJacobi) 01413 { 01414 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ; 01415 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ; 01416 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ; 01417 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ; 01418 } 01419 else 01420 { 01421 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ; 01422 } 01423 alphaBeta = alpha + beta ; 01424 for (k=1;k<=maxNumber;k++) 01425 { 01426 temp = 2.0 + alphaBeta ; 01427 nwt1 = (alpha-beta+temp*root)/2.0 ; 01428 nwt2 = 1.0 ; 01429 for (j=2;j<=nJacobi;j++) 01430 { 01431 nwt3 = nwt2 ; 01432 nwt2 = nwt1 ; 01433 temp = 2*j+alphaBeta ; 01434 a = 2*j*(j+alphaBeta)*(temp-2.0) ; 01435 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ; 01436 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ; 01437 nwt1 = (b*nwt2-c*nwt3)/a ; 01438 } 01439 nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 + 01440 2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2) / 01441 (temp*(1.0 - root*root)) ; 01442 rootTemp = root ; 01443 root = rootTemp - nwt1/nwt ; 01444 if (std::fabs(root-rootTemp) <= tolerance) 01445 { 01446 break ; 01447 } 01448 } 01449 if (k > maxNumber) 01450 { 01451 G4Exception("G4Integrator<T,F>::Jacobi(...)", "Error", 01452 FatalException, "Too many (>12) iterations."); 01453 } 01454 fAbscissa[i-1] = root ; 01455 fWeight[i-1] = 01456 std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 01457 GammaLogarithm((G4double)(beta+nJacobi)) - 01458 GammaLogarithm((G4double)(nJacobi+1.0)) - 01459 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) 01460 *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2); 01461 } 01462 01463 // 01464 // Calculation of the integral 01465 // 01466 01467 G4double integral = 0.0 ; 01468 for(i=0;i<fNumber;i++) 01469 { 01470 integral += fWeight[i]*(*f)(fAbscissa[i]) ; 01471 } 01472 delete[] fAbscissa; 01473 delete[] fWeight; 01474 return integral ; 01475 }
G4double G4Integrator< T, F >::Jacobi | ( | T * | ptrT, | |
F | f, | |||
G4double | alpha, | |||
G4double | beta, | |||
G4int | n | |||
) |
G4double G4Integrator< T, F >::Jacobi | ( | T & | typeT, | |
F | f, | |||
G4double | alpha, | |||
G4double | beta, | |||
G4int | n | |||
) |
Definition at line 1237 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().
Referenced by G4Integrator< T, F >::Jacobi().
01239 { 01240 const G4double tolerance = 1.0e-12 ; 01241 const G4double maxNumber = 12 ; 01242 G4int i, k, j ; 01243 G4double alphaBeta, alphaReduced, betaReduced, root1=0., root2=0., root3=0. ; 01244 G4double a, b, c, nwt1, nwt2, nwt3, nwt, temp, root=0., rootTemp ; 01245 01246 G4int fNumber = nJacobi ; 01247 G4double* fAbscissa = new G4double[fNumber] ; 01248 G4double* fWeight = new G4double[fNumber] ; 01249 01250 for (i=1;i<=nJacobi;i++) 01251 { 01252 if (i == 1) 01253 { 01254 alphaReduced = alpha/nJacobi ; 01255 betaReduced = beta/nJacobi ; 01256 root1 = (1.0+alpha)*(2.78002/(4.0+nJacobi*nJacobi)+ 01257 0.767999*alphaReduced/nJacobi) ; 01258 root2 = 1.0+1.48*alphaReduced+0.96002*betaReduced + 01259 0.451998*alphaReduced*alphaReduced + 01260 0.83001*alphaReduced*betaReduced ; 01261 root = 1.0-root1/root2 ; 01262 } 01263 else if (i == 2) 01264 { 01265 root1=(4.1002+alpha)/((1.0+alpha)*(1.0+0.155998*alpha)) ; 01266 root2=1.0+0.06*(nJacobi-8.0)*(1.0+0.12*alpha)/nJacobi ; 01267 root3=1.0+0.012002*beta*(1.0+0.24997*std::fabs(alpha))/nJacobi ; 01268 root -= (1.0-root)*root1*root2*root3 ; 01269 } 01270 else if (i == 3) 01271 { 01272 root1=(1.67001+0.27998*alpha)/(1.0+0.37002*alpha) ; 01273 root2=1.0+0.22*(nJacobi-8.0)/nJacobi ; 01274 root3=1.0+8.0*beta/((6.28001+beta)*nJacobi*nJacobi) ; 01275 root -= (fAbscissa[0]-root)*root1*root2*root3 ; 01276 } 01277 else if (i == nJacobi-1) 01278 { 01279 root1=(1.0+0.235002*beta)/(0.766001+0.118998*beta) ; 01280 root2=1.0/(1.0+0.639002*(nJacobi-4.0)/(1.0+0.71001*(nJacobi-4.0))) ; 01281 root3=1.0/(1.0+20.0*alpha/((7.5+alpha)*nJacobi*nJacobi)) ; 01282 root += (root-fAbscissa[nJacobi-4])*root1*root2*root3 ; 01283 } 01284 else if (i == nJacobi) 01285 { 01286 root1 = (1.0+0.37002*beta)/(1.67001+0.27998*beta) ; 01287 root2 = 1.0/(1.0+0.22*(nJacobi-8.0)/nJacobi) ; 01288 root3 = 1.0/(1.0+8.0*alpha/((6.28002+alpha)*nJacobi*nJacobi)) ; 01289 root += (root-fAbscissa[nJacobi-3])*root1*root2*root3 ; 01290 } 01291 else 01292 { 01293 root = 3.0*fAbscissa[i-2]-3.0*fAbscissa[i-3]+fAbscissa[i-4] ; 01294 } 01295 alphaBeta = alpha + beta ; 01296 for (k=1;k<=maxNumber;k++) 01297 { 01298 temp = 2.0 + alphaBeta ; 01299 nwt1 = (alpha-beta+temp*root)/2.0 ; 01300 nwt2 = 1.0 ; 01301 for (j=2;j<=nJacobi;j++) 01302 { 01303 nwt3 = nwt2 ; 01304 nwt2 = nwt1 ; 01305 temp = 2*j+alphaBeta ; 01306 a = 2*j*(j+alphaBeta)*(temp-2.0) ; 01307 b = (temp-1.0)*(alpha*alpha-beta*beta+temp*(temp-2.0)*root) ; 01308 c = 2.0*(j-1+alpha)*(j-1+beta)*temp ; 01309 nwt1 = (b*nwt2-c*nwt3)/a ; 01310 } 01311 nwt = (nJacobi*(alpha - beta - temp*root)*nwt1 + 01312 2.0*(nJacobi + alpha)*(nJacobi + beta)*nwt2)/ 01313 (temp*(1.0 - root*root)) ; 01314 rootTemp = root ; 01315 root = rootTemp - nwt1/nwt ; 01316 if (std::fabs(root-rootTemp) <= tolerance) 01317 { 01318 break ; 01319 } 01320 } 01321 if (k > maxNumber) 01322 { 01323 G4Exception("G4Integrator<T,F>::Jacobi(T,F, ...)", "Error", 01324 FatalException, "Too many (>12) iterations."); 01325 } 01326 fAbscissa[i-1] = root ; 01327 fWeight[i-1] = std::exp(GammaLogarithm((G4double)(alpha+nJacobi)) + 01328 GammaLogarithm((G4double)(beta+nJacobi)) - 01329 GammaLogarithm((G4double)(nJacobi+1.0)) - 01330 GammaLogarithm((G4double)(nJacobi + alphaBeta + 1.0))) 01331 *temp*std::pow(2.0,alphaBeta)/(nwt*nwt2) ; 01332 } 01333 01334 // 01335 // Calculation of the integral 01336 // 01337 01338 G4double integral = 0.0 ; 01339 for(i=0;i<fNumber;i++) 01340 { 01341 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ; 01342 } 01343 delete[] fAbscissa; 01344 delete[] fWeight; 01345 return integral ; 01346 }
G4double G4Integrator< T, F >::Laguerre | ( | G4double(*)(G4double) | f, | |
G4double | alpha, | |||
G4int | n | |||
) |
Definition at line 932 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().
00934 { 00935 const G4double tolerance = 1.0e-10 ; 00936 const G4int maxNumber = 12 ; 00937 G4int i, j, k ; 00938 G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ; 00939 G4double integral = 0.0 ; 00940 00941 G4int fNumber = nLaguerre ; 00942 G4double* fAbscissa = new G4double[fNumber] ; 00943 G4double* fWeight = new G4double[fNumber] ; 00944 00945 for(i=1;i<=fNumber;i++) // Loop over the desired roots 00946 { 00947 if(i == 1) 00948 { 00949 nwt = (1.0 + alpha)*(3.0 + 0.92*alpha) 00950 / (1.0 + 2.4*fNumber + 1.8*alpha) ; 00951 } 00952 else if(i == 2) 00953 { 00954 nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ; 00955 } 00956 else 00957 { 00958 cofi = i - 2 ; 00959 nwt += ((1.0+2.55*cofi)/(1.9*cofi) 00960 + 1.26*cofi*alpha/(1.0+3.5*cofi)) 00961 * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ; 00962 } 00963 for(k=1;k<=maxNumber;k++) 00964 { 00965 temp1 = 1.0 ; 00966 temp2 = 0.0 ; 00967 00968 for(j=1;j<=fNumber;j++) 00969 { 00970 temp3 = temp2 ; 00971 temp2 = temp1 ; 00972 temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ; 00973 } 00974 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ; 00975 nwt1 = nwt ; 00976 nwt = nwt1 - temp1/temp ; 00977 00978 if(std::fabs(nwt - nwt1) <= tolerance) 00979 { 00980 break ; 00981 } 00982 } 00983 if(k > maxNumber) 00984 { 00985 G4Exception("G4Integrator<T,F>::Laguerre( ...)", "Error", 00986 FatalException, "Too many (>12) iterations."); 00987 } 00988 00989 fAbscissa[i-1] = nwt ; 00990 fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 00991 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ; 00992 } 00993 00994 // 00995 // Integral evaluation 00996 // 00997 00998 for(i=0;i<fNumber;i++) 00999 { 01000 integral += fWeight[i]*(*f)(fAbscissa[i]) ; 01001 } 01002 delete[] fAbscissa; 01003 delete[] fWeight; 01004 return integral ; 01005 }
G4double G4Integrator< T, F >::Laguerre | ( | T * | ptrT, | |
F | f, | |||
G4double | alpha, | |||
G4int | n | |||
) |
G4double G4Integrator< T, F >::Laguerre | ( | T & | typeT, | |
F | f, | |||
G4double | alpha, | |||
G4int | n | |||
) |
Definition at line 840 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4Integrator< T, F >::GammaLogarithm().
Referenced by G4SynchrotronRadiationInMat::GetAngleK(), G4SynchrotronRadiationInMat::GetEnergyProbSR(), G4SynchrotronRadiationInMat::GetIntProbSR(), and G4Integrator< T, F >::Laguerre().
00842 { 00843 const G4double tolerance = 1.0e-10 ; 00844 const G4int maxNumber = 12 ; 00845 G4int i, j, k ; 00846 G4double nwt=0., nwt1, temp1, temp2, temp3, temp, cofi ; 00847 G4double integral = 0.0 ; 00848 00849 G4int fNumber = nLaguerre ; 00850 G4double* fAbscissa = new G4double[fNumber] ; 00851 G4double* fWeight = new G4double[fNumber] ; 00852 00853 for(i=1;i<=fNumber;i++) // Loop over the desired roots 00854 { 00855 if(i == 1) 00856 { 00857 nwt = (1.0 + alpha)*(3.0 + 0.92*alpha) 00858 / (1.0 + 2.4*fNumber + 1.8*alpha) ; 00859 } 00860 else if(i == 2) 00861 { 00862 nwt += (15.0 + 6.25*alpha)/(1.0 + 0.9*alpha + 2.5*fNumber) ; 00863 } 00864 else 00865 { 00866 cofi = i - 2 ; 00867 nwt += ((1.0+2.55*cofi)/(1.9*cofi) 00868 + 1.26*cofi*alpha/(1.0+3.5*cofi)) 00869 * (nwt - fAbscissa[i-3])/(1.0 + 0.3*alpha) ; 00870 } 00871 for(k=1;k<=maxNumber;k++) 00872 { 00873 temp1 = 1.0 ; 00874 temp2 = 0.0 ; 00875 00876 for(j=1;j<=fNumber;j++) 00877 { 00878 temp3 = temp2 ; 00879 temp2 = temp1 ; 00880 temp1 = ((2*j - 1 + alpha - nwt)*temp2 - (j - 1 + alpha)*temp3)/j ; 00881 } 00882 temp = (fNumber*temp1 - (fNumber +alpha)*temp2)/nwt ; 00883 nwt1 = nwt ; 00884 nwt = nwt1 - temp1/temp ; 00885 00886 if(std::fabs(nwt - nwt1) <= tolerance) 00887 { 00888 break ; 00889 } 00890 } 00891 if(k > maxNumber) 00892 { 00893 G4Exception("G4Integrator<T,F>::Laguerre(T,F, ...)", "Error", 00894 FatalException, "Too many (>12) iterations."); 00895 } 00896 00897 fAbscissa[i-1] = nwt ; 00898 fWeight[i-1] = -std::exp(GammaLogarithm(alpha + fNumber) - 00899 GammaLogarithm((G4double)fNumber))/(temp*fNumber*temp2) ; 00900 } 00901 00902 // 00903 // Integral evaluation 00904 // 00905 00906 for(i=0;i<fNumber;i++) 00907 { 00908 integral += fWeight[i]*(typeT.*f)(fAbscissa[i]) ; 00909 } 00910 delete[] fAbscissa; 00911 delete[] fWeight; 00912 return integral ; 00913 }
G4double G4Integrator< T, F >::Legendre | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 398 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4INCL::Math::pi.
00400 { 00401 G4double nwt, nwt1, temp1, temp2, temp3, temp ; 00402 G4double xDiff, xMean, dx, integral ; 00403 00404 const G4double tolerance = 1.6e-10 ; 00405 G4int i, j, k = nLegendre ; 00406 G4int fNumber = (nLegendre + 1)/2 ; 00407 00408 if(2*fNumber != k) 00409 { 00410 G4Exception("G4Integrator<T,F>::Legendre(...)", "InvalidCall", 00411 FatalException, "Invalid (odd) nLegendre in constructor."); 00412 } 00413 00414 G4double* fAbscissa = new G4double[fNumber] ; 00415 G4double* fWeight = new G4double[fNumber] ; 00416 00417 for(i=1;i<=fNumber;i++) // Loop over the desired roots 00418 { 00419 nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation 00420 00421 do // loop of Newton's method 00422 { 00423 temp1 = 1.0 ; 00424 temp2 = 0.0 ; 00425 for(j=1;j<=k;j++) 00426 { 00427 temp3 = temp2 ; 00428 temp2 = temp1 ; 00429 temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ; 00430 } 00431 temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ; 00432 nwt1 = nwt ; 00433 nwt = nwt1 - temp1/temp ; // Newton's method 00434 } 00435 while(std::fabs(nwt - nwt1) > tolerance) ; 00436 00437 fAbscissa[fNumber-i] = nwt ; 00438 fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ; 00439 } 00440 00441 // 00442 // Now we ready to get integral 00443 // 00444 00445 xMean = 0.5*(a + b) ; 00446 xDiff = 0.5*(b - a) ; 00447 integral = 0.0 ; 00448 for(i=0;i<fNumber;i++) 00449 { 00450 dx = xDiff*fAbscissa[i] ; 00451 integral += fWeight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx) ) ; 00452 } 00453 delete[] fAbscissa; 00454 delete[] fWeight; 00455 00456 return integral *= xDiff ; 00457 }
G4double G4Integrator< T, F >::Legendre | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 387 of file G4Integrator.icc.
References G4Integrator< T, F >::Legendre().
00389 { 00390 return Legendre(*ptrT,f,a,b,nLegendre) ; 00391 }
G4double G4Integrator< T, F >::Legendre | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 321 of file G4Integrator.icc.
References FatalException, G4Exception(), and G4INCL::Math::pi.
Referenced by G4Integrator< T, F >::Legendre().
00323 { 00324 G4double nwt, nwt1, temp1, temp2, temp3, temp ; 00325 G4double xDiff, xMean, dx, integral ; 00326 00327 const G4double tolerance = 1.6e-10 ; 00328 G4int i, j, k = nLegendre ; 00329 G4int fNumber = (nLegendre + 1)/2 ; 00330 00331 if(2*fNumber != k) 00332 { 00333 G4Exception("G4Integrator<T,F>::Legendre(T&,F, ...)", "InvalidCall", 00334 FatalException, "Invalid (odd) nLegendre in constructor."); 00335 } 00336 00337 G4double* fAbscissa = new G4double[fNumber] ; 00338 G4double* fWeight = new G4double[fNumber] ; 00339 00340 for(i=1;i<=fNumber;i++) // Loop over the desired roots 00341 { 00342 nwt = std::cos(CLHEP::pi*(i - 0.25)/(k + 0.5)) ; // Initial root approximation 00343 00344 do // loop of Newton's method 00345 { 00346 temp1 = 1.0 ; 00347 temp2 = 0.0 ; 00348 for(j=1;j<=k;j++) 00349 { 00350 temp3 = temp2 ; 00351 temp2 = temp1 ; 00352 temp1 = ((2.0*j - 1.0)*nwt*temp2 - (j - 1.0)*temp3)/j ; 00353 } 00354 temp = k*(nwt*temp1 - temp2)/(nwt*nwt - 1.0) ; 00355 nwt1 = nwt ; 00356 nwt = nwt1 - temp1/temp ; // Newton's method 00357 } 00358 while(std::fabs(nwt - nwt1) > tolerance) ; 00359 00360 fAbscissa[fNumber-i] = nwt ; 00361 fWeight[fNumber-i] = 2.0/((1.0 - nwt*nwt)*temp*temp) ; 00362 } 00363 00364 // 00365 // Now we ready to get integral 00366 // 00367 00368 xMean = 0.5*(a + b) ; 00369 xDiff = 0.5*(b - a) ; 00370 integral = 0.0 ; 00371 for(i=0;i<fNumber;i++) 00372 { 00373 dx = xDiff*fAbscissa[i] ; 00374 integral += fWeight[i]*( (typeT.*f)(xMean + dx) + 00375 (typeT.*f)(xMean - dx) ) ; 00376 } 00377 delete[] fAbscissa; 00378 delete[] fWeight; 00379 return integral *= xDiff ; 00380 }
G4double G4Integrator< T, F >::Legendre10 | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b | |||
) |
Definition at line 509 of file G4Integrator.icc.
00511 { 00512 G4int i ; 00513 G4double xDiff, xMean, dx, integral ; 00514 00515 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916 00516 00517 static G4double abscissa[] = { 0.148874338981631, 0.433395394129247, 00518 0.679409568299024, 0.865063366688985, 00519 0.973906528517172 } ; 00520 00521 static G4double weight[] = { 0.295524224714753, 0.269266719309996, 00522 0.219086362515982, 0.149451349150581, 00523 0.066671344308688 } ; 00524 xMean = 0.5*(a + b) ; 00525 xDiff = 0.5*(b - a) ; 00526 integral = 0.0 ; 00527 for(i=0;i<5;i++) 00528 { 00529 dx = xDiff*abscissa[i] ; 00530 integral += weight[i]*( (*f)(xMean + dx) + (*f)(xMean - dx)) ; 00531 } 00532 return integral *= xDiff ; 00533 }
G4double G4Integrator< T, F >::Legendre10 | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) |
Definition at line 499 of file G4Integrator.icc.
References G4Integrator< T, F >::Legendre10().
00500 { 00501 return Legendre10(*ptrT,f,a,b) ; 00502 }
G4double G4Integrator< T, F >::Legendre10 | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) |
Definition at line 469 of file G4Integrator.icc.
Referenced by G4VXTRenergyLoss::BuildAngleForEnergyBank(), G4NuclNuclDiffuseElastic::BuildAngleTable(), G4DiffuseElastic::BuildAngleTable(), G4VXTRenergyLoss::BuildEnergyTable(), G4InitXscPAI::IntegralCherenkov(), G4InitXscPAI::IntegralPAIdEdx(), G4InitXscPAI::IntegralPAIxSection(), G4InitXscPAI::IntegralPlasmon(), G4Integrator< T, F >::Legendre10(), G4NuclNuclDiffuseElastic::SampleThetaCMS(), G4DiffuseElastic::SampleThetaCMS(), G4NuclNuclDiffuseElastic::TestAngleTable(), G4DiffuseElastic::TestAngleTable(), and G4VXTRenergyLoss::XTRNSpectralDensity().
00470 { 00471 G4int i ; 00472 G4double xDiff, xMean, dx, integral ; 00473 00474 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 916 00475 00476 static G4double abscissa[] = { 0.148874338981631, 0.433395394129247, 00477 0.679409568299024, 0.865063366688985, 00478 0.973906528517172 } ; 00479 00480 static G4double weight[] = { 0.295524224714753, 0.269266719309996, 00481 0.219086362515982, 0.149451349150581, 00482 0.066671344308688 } ; 00483 xMean = 0.5*(a + b) ; 00484 xDiff = 0.5*(b - a) ; 00485 integral = 0.0 ; 00486 for(i=0;i<5;i++) 00487 { 00488 dx = xDiff*abscissa[i] ; 00489 integral += weight[i]*( (typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ; 00490 } 00491 return integral *= xDiff ; 00492 }
G4double G4Integrator< T, F >::Legendre96 | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b | |||
) |
Definition at line 647 of file G4Integrator.icc.
00649 { 00650 G4int i ; 00651 G4double xDiff, xMean, dx, integral ; 00652 00653 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919 00654 00655 static G4double 00656 abscissa[] = { 00657 0.016276744849602969579, 0.048812985136049731112, 00658 0.081297495464425558994, 0.113695850110665920911, 00659 0.145973714654896941989, 0.178096882367618602759, // 6 00660 00661 0.210031310460567203603, 0.241743156163840012328, 00662 0.273198812591049141487, 0.304364944354496353024, 00663 0.335208522892625422616, 0.365696861472313635031, // 12 00664 00665 0.395797649828908603285, 0.425478988407300545365, 00666 0.454709422167743008636, 0.483457973920596359768, 00667 0.511694177154667673586, 0.539388108324357436227, // 18 00668 00669 0.566510418561397168404, 0.593032364777572080684, 00670 0.618925840125468570386, 0.644163403784967106798, 00671 0.668718310043916153953, 0.692564536642171561344, // 24 00672 00673 0.715676812348967626225, 0.738030643744400132851, 00674 0.759602341176647498703, 0.780369043867433217604, 00675 0.800308744139140817229, 0.819400310737931675539, // 30 00676 00677 0.837623511228187121494, 0.854959033434601455463, 00678 0.871388505909296502874, 0.886894517402420416057, 00679 0.901460635315852341319, 0.915071423120898074206, // 36 00680 00681 0.927712456722308690965, 0.939370339752755216932, 00682 0.950032717784437635756, 0.959688291448742539300, 00683 0.968326828463264212174, 0.975939174585136466453, // 42 00684 00685 0.982517263563014677447, 0.988054126329623799481, 00686 0.992543900323762624572, 0.995981842987209290650, 00687 0.998364375863181677724, 0.999689503883230766828 // 48 00688 } ; 00689 00690 static G4double 00691 weight[] = { 00692 0.032550614492363166242, 0.032516118713868835987, 00693 0.032447163714064269364, 0.032343822568575928429, 00694 0.032206204794030250669, 0.032034456231992663218, // 6 00695 00696 0.031828758894411006535, 0.031589330770727168558, 00697 0.031316425596862355813, 0.031010332586313837423, 00698 0.030671376123669149014, 0.030299915420827593794, // 12 00699 00700 0.029896344136328385984, 0.029461089958167905970, 00701 0.028994614150555236543, 0.028497411065085385646, 00702 0.027970007616848334440, 0.027412962726029242823, // 18 00703 00704 0.026826866725591762198, 0.026212340735672413913, 00705 0.025570036005349361499, 0.024900633222483610288, 00706 0.024204841792364691282, 0.023483399085926219842, // 24 00707 00708 0.022737069658329374001, 0.021966644438744349195, 00709 0.021172939892191298988, 0.020356797154333324595, 00710 0.019519081140145022410, 0.018660679627411467385, // 30 00711 00712 0.017782502316045260838, 0.016885479864245172450, 00713 0.015970562902562291381, 0.015038721026994938006, 00714 0.014090941772314860916, 0.013128229566961572637, // 36 00715 00716 0.012151604671088319635, 0.011162102099838498591, 00717 0.010160770535008415758, 0.009148671230783386633, 00718 0.008126876925698759217, 0.007096470791153865269, // 42 00719 00720 0.006058545504235961683, 0.005014202742927517693, 00721 0.003964554338444686674, 0.002910731817934946408, 00722 0.001853960788946921732, 0.000796792065552012429 // 48 00723 } ; 00724 xMean = 0.5*(a + b) ; 00725 xDiff = 0.5*(b - a) ; 00726 integral = 0.0 ; 00727 for(i=0;i<48;i++) 00728 { 00729 dx = xDiff*abscissa[i] ; 00730 integral += weight[i]*((*f)(xMean + dx) + (*f)(xMean - dx)) ; 00731 } 00732 return integral *= xDiff ; 00733 }
G4double G4Integrator< T, F >::Legendre96 | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) |
Definition at line 637 of file G4Integrator.icc.
References G4Integrator< T, F >::Legendre96().
00638 { 00639 return Legendre96(*ptrT,f,a,b) ; 00640 }
G4double G4Integrator< T, F >::Legendre96 | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b | |||
) |
Definition at line 545 of file G4Integrator.icc.
Referenced by G4VXTRenergyLoss::BuildGlobalAngleTable(), G4NuclNuclDiffuseElastic::GetCint(), G4NuclNuclDiffuseElastic::GetErfInt(), G4NuclNuclDiffuseElastic::GetSint(), G4NuclNuclDiffuseElastic::IntegralElasticProb(), G4DiffuseElastic::IntegralElasticProb(), G4Integrator< T, F >::Legendre96(), G4NuclNuclDiffuseElastic::SampleThetaCMS(), G4DiffuseElastic::SampleThetaCMS(), G4VXTRenergyLoss::SpectralXTRdEdx(), G4NuclNuclDiffuseElastic::TestAngleTable(), G4DiffuseElastic::TestAngleTable(), G4VXTRenergyLoss::XTRNAngleDensity(), and G4VXTRenergyLoss::XTRNSpectralDensity().
00546 { 00547 G4int i ; 00548 G4double xDiff, xMean, dx, integral ; 00549 00550 // From Abramowitz M., Stegan I.A. 1964 , Handbook of Math... , p. 919 00551 00552 static G4double 00553 abscissa[] = { 00554 0.016276744849602969579, 0.048812985136049731112, 00555 0.081297495464425558994, 0.113695850110665920911, 00556 0.145973714654896941989, 0.178096882367618602759, // 6 00557 00558 0.210031310460567203603, 0.241743156163840012328, 00559 0.273198812591049141487, 0.304364944354496353024, 00560 0.335208522892625422616, 0.365696861472313635031, // 12 00561 00562 0.395797649828908603285, 0.425478988407300545365, 00563 0.454709422167743008636, 0.483457973920596359768, 00564 0.511694177154667673586, 0.539388108324357436227, // 18 00565 00566 0.566510418561397168404, 0.593032364777572080684, 00567 0.618925840125468570386, 0.644163403784967106798, 00568 0.668718310043916153953, 0.692564536642171561344, // 24 00569 00570 0.715676812348967626225, 0.738030643744400132851, 00571 0.759602341176647498703, 0.780369043867433217604, 00572 0.800308744139140817229, 0.819400310737931675539, // 30 00573 00574 0.837623511228187121494, 0.854959033434601455463, 00575 0.871388505909296502874, 0.886894517402420416057, 00576 0.901460635315852341319, 0.915071423120898074206, // 36 00577 00578 0.927712456722308690965, 0.939370339752755216932, 00579 0.950032717784437635756, 0.959688291448742539300, 00580 0.968326828463264212174, 0.975939174585136466453, // 42 00581 00582 0.982517263563014677447, 0.988054126329623799481, 00583 0.992543900323762624572, 0.995981842987209290650, 00584 0.998364375863181677724, 0.999689503883230766828 // 48 00585 } ; 00586 00587 static G4double 00588 weight[] = { 00589 0.032550614492363166242, 0.032516118713868835987, 00590 0.032447163714064269364, 0.032343822568575928429, 00591 0.032206204794030250669, 0.032034456231992663218, // 6 00592 00593 0.031828758894411006535, 0.031589330770727168558, 00594 0.031316425596862355813, 0.031010332586313837423, 00595 0.030671376123669149014, 0.030299915420827593794, // 12 00596 00597 0.029896344136328385984, 0.029461089958167905970, 00598 0.028994614150555236543, 0.028497411065085385646, 00599 0.027970007616848334440, 0.027412962726029242823, // 18 00600 00601 0.026826866725591762198, 0.026212340735672413913, 00602 0.025570036005349361499, 0.024900633222483610288, 00603 0.024204841792364691282, 0.023483399085926219842, // 24 00604 00605 0.022737069658329374001, 0.021966644438744349195, 00606 0.021172939892191298988, 0.020356797154333324595, 00607 0.019519081140145022410, 0.018660679627411467385, // 30 00608 00609 0.017782502316045260838, 0.016885479864245172450, 00610 0.015970562902562291381, 0.015038721026994938006, 00611 0.014090941772314860916, 0.013128229566961572637, // 36 00612 00613 0.012151604671088319635, 0.011162102099838498591, 00614 0.010160770535008415758, 0.009148671230783386633, 00615 0.008126876925698759217, 0.007096470791153865269, // 42 00616 00617 0.006058545504235961683, 0.005014202742927517693, 00618 0.003964554338444686674, 0.002910731817934946408, 00619 0.001853960788946921732, 0.000796792065552012429 // 48 00620 } ; 00621 xMean = 0.5*(a + b) ; 00622 xDiff = 0.5*(b - a) ; 00623 integral = 0.0 ; 00624 for(i=0;i<48;i++) 00625 { 00626 dx = xDiff*abscissa[i] ; 00627 integral += weight[i]*((typeT.*f)(xMean + dx) + (typeT.*f)(xMean - dx)) ; 00628 } 00629 return integral *= xDiff ; 00630 }
G4double G4Integrator< T, F >::Simpson | ( | G4double(*)(G4double) | f, | |
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 105 of file G4Integrator.icc.
00109 { 00110 G4int i ; 00111 G4double step = (xFinal - xInitial)/iterationNumber ; 00112 G4double x = xInitial ; 00113 G4double xPlus = xInitial + 0.5*step ; 00114 G4double mean = ( (*f)(xInitial) + (*f)(xFinal) )*0.5 ; 00115 G4double sum = (*f)(xPlus) ; 00116 00117 for(i=1;i<iterationNumber;i++) 00118 { 00119 x += step ; 00120 xPlus += step ; 00121 mean += (*f)(x) ; 00122 sum += (*f)(xPlus) ; 00123 } 00124 mean += 2.0*sum ; 00125 00126 return mean*step/3.0 ; 00127 }
G4double G4Integrator< T, F >::Simpson | ( | T * | ptrT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 73 of file G4Integrator.icc.
00078 { 00079 G4int i ; 00080 G4double step = (xFinal - xInitial)/iterationNumber ; 00081 G4double x = xInitial ; 00082 G4double xPlus = xInitial + 0.5*step ; 00083 G4double mean = ( (ptrT->*f)(xInitial) + (ptrT->*f)(xFinal) )*0.5 ; 00084 G4double sum = (ptrT->*f)(xPlus) ; 00085 00086 for(i=1;i<iterationNumber;i++) 00087 { 00088 x += step ; 00089 xPlus += step ; 00090 mean += (ptrT->*f)(x) ; 00091 sum += (ptrT->*f)(xPlus) ; 00092 } 00093 mean += 2.0*sum ; 00094 00095 return mean*step/3.0 ; 00096 }
G4double G4Integrator< T, F >::Simpson | ( | T & | typeT, | |
F | f, | |||
G4double | a, | |||
G4double | b, | |||
G4int | n | |||
) |
Definition at line 42 of file G4Integrator.icc.
Referenced by G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(), G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(), G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(), and G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond().
00047 { 00048 G4int i ; 00049 G4double step = (xFinal - xInitial)/iterationNumber ; 00050 G4double x = xInitial ; 00051 G4double xPlus = xInitial + 0.5*step ; 00052 G4double mean = ( (typeT.*f)(xInitial) + (typeT.*f)(xFinal) )*0.5 ; 00053 G4double sum = (typeT.*f)(xPlus) ; 00054 00055 for(i=1;i<iterationNumber;i++) 00056 { 00057 x += step ; 00058 xPlus += step ; 00059 mean += (typeT.*f)(x) ; 00060 sum += (typeT.*f)(xPlus) ; 00061 } 00062 mean += 2.0*sum ; 00063 00064 return mean*step/3.0 ; 00065 }