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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00060
00061 #include "G4Bessel.hh"
00062 #include "G4PhysicalConstants.hh"
00064
00065 G4Bessel::G4Bessel ()
00066 {;}
00068
00069 G4Bessel::~G4Bessel ()
00070 {;}
00072
00073 G4double G4Bessel::I0 (G4double x)
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 }
00108
00109 G4double G4Bessel::K0 (G4double x)
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 }
00142
00143 G4double G4Bessel::I1 (G4double x)
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 }
00180
00181 G4double G4Bessel::K1 (G4double x)
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 }
00214
00215 G4double G4Bessel::pI0 (G4double x)
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 }
00257
00258 G4double G4Bessel::pI1 (G4double x)
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 }
00302
00303 G4double G4Bessel::pK0 (G4double x)
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 }
00343
00344 G4double G4Bessel::pK1 (G4double x)
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 }
00354