G4Bessel Class Reference

#include <G4Bessel.hh>


Public Member Functions

 G4Bessel ()
 ~G4Bessel ()
G4double I0 (G4double)
G4double I1 (G4double)
G4double K0 (G4double)
G4double K1 (G4double)
G4double pI0 (G4double)
G4double pI1 (G4double)
G4double pK0 (G4double)
G4double pK1 (G4double)

Detailed Description

Definition at line 66 of file G4Bessel.hh.

Constructor & Destructor Documentation

G4Bessel::G4Bessel (  ) 

Definition at line 65 of file G4Bessel.cc.

00066 {;}

G4Bessel::~G4Bessel (  ) 

Definition at line 69 of file G4Bessel.cc.

00070 {;}

Member Function Documentation

G4double G4Bessel::I0 ( G4double   ) 

Definition at line 73 of file G4Bessel.cc.

Referenced by K0().

00074 {
00075   const G4double P1 = 1.0,
00076                  P2 = 3.5156229,
00077                  P3 = 3.0899424,
00078                  P4 = 1.2067492,
00079                  P5 = 0.2659732,
00080                  P6 = 0.0360768,
00081                  P7 = 0.0045813;
00082   const G4double Q1 = 0.39894228,
00083                  Q2 = 0.01328592,
00084                  Q3 = 0.00225319,
00085                  Q4 =-0.00157565,
00086                  Q5 = 0.00916281,
00087                  Q6 =-0.02057706,
00088                  Q7 = 0.02635537,
00089                  Q8 =-0.01647633,
00090                  Q9 = 0.00392377;
00091   
00092   G4double I = 0.0;
00093   if (std::fabs(x) < 3.75)
00094   {
00095     G4double y = std::pow(x/3.75, 2.0);
00096     I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
00097   }
00098   else
00099   {
00100     G4double ax = std::fabs(x);
00101     G4double y  = 3.75/ax;
00102     I  = std::exp(ax) / std::sqrt(ax) *
00103       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
00104   }
00105   return I;
00106 }

G4double G4Bessel::I1 ( G4double   ) 

Definition at line 143 of file G4Bessel.cc.

Referenced by K1().

00144 {
00145   const G4double P1 = 0.5,
00146                  P2 = 0.87890594,
00147                  P3 = 0.51498869,
00148                  P4 = 0.15084934,
00149                  P5 = 0.02658733,
00150                  P6 = 0.00301532,
00151                  P7 = 0.00032411;
00152   const G4double Q1 = 0.39894228,
00153                  Q2 =-0.03988024,
00154                  Q3 =-0.00362018,
00155                  Q4 = 0.00163801,
00156                  Q5 =-0.01031555,
00157                  Q6 = 0.02282967,
00158                  Q7 =-0.02895312,
00159                  Q8 = 0.01787654,
00160                  Q9 =-0.00420059;
00161 
00162   G4double I = 0.0;
00163   if (std::fabs(x) < 3.75)
00164   {
00165     G4double ax = std::fabs(x);
00166     G4double y = std::pow(x/3.75, 2.0);
00167     I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
00168   }
00169   else
00170   {
00171     G4double ax = std::fabs(x);
00172     G4double y  = 3.75/ax;
00173     I  = std::exp(ax) / std::sqrt(ax) *
00174       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
00175   }
00176   if (x < 0.0) I = -I;
00177   return I;
00178 }

G4double G4Bessel::K0 ( G4double   ) 

Definition at line 109 of file G4Bessel.cc.

References I0().

00110 {
00111   const G4double P1 =-0.57721566,
00112                  P2 = 0.42278420,
00113                  P3 = 0.23069756,
00114                  P4 = 0.03488590,
00115                  P5 = 0.00262698,
00116                  P6 = 0.00010750,
00117                  P7 = 0.00000740;
00118   const G4double Q1 = 1.25331414,
00119                  Q2 =-0.07832358,
00120                  Q3 = 0.02189568,
00121                  Q4 =-0.01062446,
00122                  Q5 = 0.00587872,
00123                  Q6 =-0.00251540,
00124                  Q7 = 0.00053208;
00125 
00126   G4double K = 0.0;
00127   if (x <= 2.0)
00128   {
00129     G4double y = x * x / 4.0;
00130     K = (-std::log(x/2.0)) * I0(x) +
00131       P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
00132   }
00133   else
00134   {
00135     G4double y = 2.0 / x;
00136     K = std::exp(-x)  / std::sqrt(x) *
00137       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
00138   }
00139   return K;
00140 }

G4double G4Bessel::K1 ( G4double   ) 

Definition at line 181 of file G4Bessel.cc.

References I1().

00182 {
00183   const G4double P1 = 1.0,
00184                  P2 = 0.15443144,
00185                  P3 =-0.67278579,
00186                  P4 =-0.18156897,
00187                  P5 =-0.01919402,
00188                  P6 =-0.00110404,
00189                  P7 =-0.00004686;
00190   const G4double Q1 = 1.25331414,
00191                  Q2 = 0.23498619,
00192                  Q3 =-0.03655620,
00193                  Q4 = 0.01504268,
00194                  Q5 =-0.00780353,
00195                  Q6 = 0.00325614,
00196                  Q7 =-0.00068245;
00197 
00198   G4double K = 0.0;
00199   if (x <= 2.0)
00200   {
00201     G4double y = x * x / 4.0;
00202     K = std::log(x/2.0)*I1(x) + 1.0/x *
00203       (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
00204   }
00205   else
00206   {
00207     G4double y = 2.0 / x;
00208     K = std::exp(-x) / std::sqrt(x) *
00209       (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
00210   }
00211   return K;
00212 }

G4double G4Bessel::pI0 ( G4double   ) 

Definition at line 215 of file G4Bessel.cc.

References A10, and A11.

Referenced by pK0(), and pK1().

00216 {
00217   const G4double A0  = 0.1250000000000E+00,
00218                  A1  = 7.0312500000000E-02,
00219                  A2  = 7.3242187500000E-02,
00220                  A3  = 1.1215209960938E-01,
00221                  A4  = 2.2710800170898E-01,
00222                  A5  = 5.7250142097473E-01,
00223                  A6  = 1.7277275025845E+00,
00224                  A7  = 6.0740420012735E+00,
00225                  A8  = 2.4380529699556E+01,
00226                  A9  = 1.1001714026925E+02,
00227                  A10 = 5.5133589612202E+02,
00228                  A11 = 3.0380905109224E+03;
00229 
00230   G4double I = 0.0;
00231   if (x == 0.0)
00232   {
00233     I = 1.0;
00234   }
00235   else if (x < 18.0)
00236   {
00237     I          = 1.0;
00238     G4double y = x * x;
00239     G4double q = 1.0;
00240     for (G4int i=1; i<101; i++)
00241     {
00242       q *= 0.25 * y / i / i;
00243       I += q;
00244       if (std::abs(q/I) < 1.0E-15) break;
00245     }
00246   }
00247   else
00248   {
00249     G4double y = 1.0 / x;
00250     I = std::exp(x) / std::sqrt(twopi*x) *
00251       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
00252   }
00253 
00254   return I;
00255 }

G4double G4Bessel::pI1 ( G4double   ) 

Definition at line 258 of file G4Bessel.cc.

References A10, and A11.

Referenced by pK1().

00259 {
00260   const G4double A0  = -0.3750000000000E+00,
00261                  A1  = -1.1718750000000E-01,
00262                  A2  = -1.0253906250000E-01,
00263                  A3  = -1.4419555664063E-01,
00264                  A4  = -2.775764465332E-01,
00265                  A5  = -6.7659258842468E-01,
00266                  A6  = -1.9935317337513E+00,
00267                  A7  = -6.8839142681099E+00,
00268                  A8  = -2.7248827311269E+01,
00269                  A9  = -1.2159789187654E+02,
00270                  A10 = -6.0384407670507E+02,
00271                  A11 = -3.3022722944809E+03;
00272 
00273   G4double I = 0.0;
00274   if (x == 0.0)
00275   {
00276     I = 0.0;
00277   }
00278   else if (x < 18.0)
00279   {
00280     I          = 1.0;
00281     G4double y = x * x;
00282     G4double q = 1.0;
00283     for (G4int i=1; i<101; i++)
00284     {
00285       q *= 0.25 * y / i / (i+1.0);
00286       I += q;
00287       if (std::abs(q/I) < 1.0E-15) break;
00288     }
00289     I *= 0.5 * x;
00290     
00291   }
00292   else
00293   {
00294     G4double y = 1.0 / x;
00295     I = std::exp(x) / std::sqrt(twopi*x) *
00296       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11))))))))))));
00297   }
00298 
00299   return I;
00300 }

G4double G4Bessel::pK0 ( G4double   ) 

Definition at line 303 of file G4Bessel.cc.

References pI0().

Referenced by pK1().

00304 {
00305   const G4double A0 = 0.1250000000000E+00,
00306                  A1 = 0.2109375000000E+00,
00307                  A2 = 1.0986328125000E+00,
00308                  A3 = 1.1775970458984E+01,
00309                  A4 = 2.1461706161499E+02,
00310                  A5 = 5.9511522710323E+03,
00311                  A6 = 2.3347645606175E+05,
00312                  A7 = 1.2312234987631E+07;
00313 
00314   G4double K = 0.0;
00315   if (x == 0.0)
00316   {
00317     K = 1.0E+307;
00318   }
00319   else if (x < 9.0)
00320   {
00321     G4double y = x * x;
00322     G4double C = -std::log(x/2.0) - 0.5772156649015329;
00323     G4double q = 1.0;
00324     G4double t = 0.0;
00325     for (G4int i=1; i<51; i++)
00326     {
00327       q *= 0.25 * y / i / i;
00328       t += 1.0 / i ;
00329       K += q * (t+C);
00330     }
00331     K += C;
00332   }
00333   else
00334   {
00335     G4double y = 1.0 / x / x;
00336     K = 0.5 / x / pI0(x) *
00337       (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7))))))));
00338   }
00339   
00340   return K;
00341 }

G4double G4Bessel::pK1 ( G4double   ) 

Definition at line 344 of file G4Bessel.cc.

References pI0(), pI1(), and pK0().

00345 {
00346   G4double K = 0.0;
00347   if (x == 0.0)
00348     K = 1.0E+307;
00349   else
00350     K = (1.0/x - pI1(x)*pK0(x)) / pI0(x);
00351   return K;
00352 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:31 2013 for Geant4 by  doxygen 1.4.7