Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4Bessel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * *
21// * Parts of this code which have been developed by QinetiQ Ltd *
22// * under contract to the European Space Agency (ESA) are the *
23// * intellectual property of ESA. Rights to use, copy, modify and *
24// * redistribute this software for general public use are granted *
25// * in compliance with any licensing, distribution and development *
26// * policy adopted by the Geant4 Collaboration. This code has been *
27// * written by QinetiQ Ltd for the European Space Agency, under ESA *
28// * contract 17191/03/NL/LvH (Aurora Programme). *
29// * *
30// * By using, copying, modifying or distributing the software (or *
31// * any work based on the software) you agree to acknowledge its *
32// * use in resulting scientific publications, and indicate your *
33// * acceptance of all terms of the Geant4 Software license. *
34// ********************************************************************
35//
36// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37//
38// MODULE: G4Bessel.cc
39//
40// Version: B.1
41// Date: 15/04/04
42// Author: P R Truscott
43// Organisation: QinetiQ Ltd, UK
44// Customer: ESA/ESTEC, NOORDWIJK
45// Contract: 17191/03/NL/LvH
46//
47// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48//
49// CHANGE HISTORY
50// --------------
51//
52// 18 Noevmber 2003, P R Truscott, QinetiQ Ltd, UK
53// Created.
54//
55// 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56// Beta release
57//
58// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59////////////////////////////////////////////////////////////////////////////////
60//
61#include "G4Bessel.hh"
62#include "G4PhysicalConstants.hh"
63////////////////////////////////////////////////////////////////////////////////
64//
65G4Bessel::G4Bessel ()
66{;}
67////////////////////////////////////////////////////////////////////////////////
68//
69G4Bessel::~G4Bessel ()
70{;}
71////////////////////////////////////////////////////////////////////////////////
72//
73G4double G4Bessel::I0 (G4double x)
74{
75 const G4double P1 = 1.0,
76 P2 = 3.5156229,
77 P3 = 3.0899424,
78 P4 = 1.2067492,
79 P5 = 0.2659732,
80 P6 = 0.0360768,
81 P7 = 0.0045813;
82 const G4double Q1 = 0.39894228,
83 Q2 = 0.01328592,
84 Q3 = 0.00225319,
85 Q4 =-0.00157565,
86 Q5 = 0.00916281,
87 Q6 =-0.02057706,
88 Q7 = 0.02635537,
89 Q8 =-0.01647633,
90 Q9 = 0.00392377;
91
92 G4double I = 0.0;
93 if (std::fabs(x) < 3.75)
94 {
95 G4double y = std::pow(x/3.75, 2.0);
96 I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
97 }
98 else
99 {
100 G4double ax = std::fabs(x);
101 G4double y = 3.75/ax;
102 I = std::exp(ax) / std::sqrt(ax) *
103 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
104 }
105 return I;
106}
107////////////////////////////////////////////////////////////////////////////////
108//
109G4double G4Bessel::K0 (G4double x)
110{
111 const G4double P1 =-0.57721566,
112 P2 = 0.42278420,
113 P3 = 0.23069756,
114 P4 = 0.03488590,
115 P5 = 0.00262698,
116 P6 = 0.00010750,
117 P7 = 0.00000740;
118 const G4double Q1 = 1.25331414,
119 Q2 =-0.07832358,
120 Q3 = 0.02189568,
121 Q4 =-0.01062446,
122 Q5 = 0.00587872,
123 Q6 =-0.00251540,
124 Q7 = 0.00053208;
125
126 G4double K = 0.0;
127 if (x <= 2.0)
128 {
129 G4double y = x * x / 4.0;
130 K = (-std::log(x/2.0)) * I0(x) +
131 P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
132 }
133 else
134 {
135 G4double y = 2.0 / x;
136 K = std::exp(-x) / std::sqrt(x) *
137 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
138 }
139 return K;
140}
141////////////////////////////////////////////////////////////////////////////////
142//
143G4double G4Bessel::I1 (G4double x)
144{
145 const G4double P1 = 0.5,
146 P2 = 0.87890594,
147 P3 = 0.51498869,
148 P4 = 0.15084934,
149 P5 = 0.02658733,
150 P6 = 0.00301532,
151 P7 = 0.00032411;
152 const G4double Q1 = 0.39894228,
153 Q2 =-0.03988024,
154 Q3 =-0.00362018,
155 Q4 = 0.00163801,
156 Q5 =-0.01031555,
157 Q6 = 0.02282967,
158 Q7 =-0.02895312,
159 Q8 = 0.01787654,
160 Q9 =-0.00420059;
161
162 G4double I = 0.0;
163 if (std::fabs(x) < 3.75)
164 {
165 G4double ax = std::fabs(x);
166 G4double y = std::pow(x/3.75, 2.0);
167 I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
168 }
169 else
170 {
171 G4double ax = std::fabs(x);
172 G4double y = 3.75/ax;
173 I = std::exp(ax) / std::sqrt(ax) *
174 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
175 }
176 if (x < 0.0) I = -I;
177 return I;
178}
179////////////////////////////////////////////////////////////////////////////////
180//
181G4double G4Bessel::K1 (G4double x)
182{
183 const G4double P1 = 1.0,
184 P2 = 0.15443144,
185 P3 =-0.67278579,
186 P4 =-0.18156897,
187 P5 =-0.01919402,
188 P6 =-0.00110404,
189 P7 =-0.00004686;
190 const G4double Q1 = 1.25331414,
191 Q2 = 0.23498619,
192 Q3 =-0.03655620,
193 Q4 = 0.01504268,
194 Q5 =-0.00780353,
195 Q6 = 0.00325614,
196 Q7 =-0.00068245;
197
198 G4double K = 0.0;
199 if (x <= 2.0)
200 {
201 G4double y = x * x / 4.0;
202 K = std::log(x/2.0)*I1(x) + 1.0/x *
203 (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
204 }
205 else
206 {
207 G4double y = 2.0 / x;
208 K = std::exp(-x) / std::sqrt(x) *
209 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
210 }
211 return K;
212}
213////////////////////////////////////////////////////////////////////////////////
214//
215G4double G4Bessel::pI0 (G4double x)
216{
217 const G4double A0 = 0.1250000000000E+00,
218 A1 = 7.0312500000000E-02,
219 A2 = 7.3242187500000E-02,
220 A3 = 1.1215209960938E-01,
221 A4 = 2.2710800170898E-01,
222 A5 = 5.7250142097473E-01,
223 A6 = 1.7277275025845E+00,
224 A7 = 6.0740420012735E+00,
225 A8 = 2.4380529699556E+01,
226 A9 = 1.1001714026925E+02,
227 A10 = 5.5133589612202E+02,
228 A11 = 3.0380905109224E+03;
229
230 G4double I = 0.0;
231 if (x == 0.0)
232 {
233 I = 1.0;
234 }
235 else if (x < 18.0)
236 {
237 I = 1.0;
238 G4double y = x * x;
239 G4double q = 1.0;
240 for (G4int i=1; i<101; i++)
241 {
242 q *= 0.25 * y / i / i;
243 I += q;
244 if (std::abs(q/I) < 1.0E-15) break;
245 }
246 }
247 else
248 {
249 G4double y = 1.0 / x;
250 I = std::exp(x) / std::sqrt(twopi*x) *
251 (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))))))))))));
252 }
253
254 return I;
255}
256////////////////////////////////////////////////////////////////////////////////
257//
258G4double G4Bessel::pI1 (G4double x)
259{
260 const G4double A0 = -0.3750000000000E+00,
261 A1 = -1.1718750000000E-01,
262 A2 = -1.0253906250000E-01,
263 A3 = -1.4419555664063E-01,
264 A4 = -2.775764465332E-01,
265 A5 = -6.7659258842468E-01,
266 A6 = -1.9935317337513E+00,
267 A7 = -6.8839142681099E+00,
268 A8 = -2.7248827311269E+01,
269 A9 = -1.2159789187654E+02,
270 A10 = -6.0384407670507E+02,
271 A11 = -3.3022722944809E+03;
272
273 G4double I = 0.0;
274 if (x == 0.0)
275 {
276 I = 0.0;
277 }
278 else if (x < 18.0)
279 {
280 I = 1.0;
281 G4double y = x * x;
282 G4double q = 1.0;
283 for (G4int i=1; i<101; i++)
284 {
285 q *= 0.25 * y / i / (i+1.0);
286 I += q;
287 if (std::abs(q/I) < 1.0E-15) break;
288 }
289 I *= 0.5 * x;
290
291 }
292 else
293 {
294 G4double y = 1.0 / x;
295 I = std::exp(x) / std::sqrt(twopi*x) *
296 (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))))))))))));
297 }
298
299 return I;
300}
301////////////////////////////////////////////////////////////////////////////////
302//
303G4double G4Bessel::pK0 (G4double x)
304{
305 const G4double A0 = 0.1250000000000E+00,
306 A1 = 0.2109375000000E+00,
307 A2 = 1.0986328125000E+00,
308 A3 = 1.1775970458984E+01,
309 A4 = 2.1461706161499E+02,
310 A5 = 5.9511522710323E+03,
311 A6 = 2.3347645606175E+05,
312 A7 = 1.2312234987631E+07;
313
314 G4double K = 0.0;
315 if (x == 0.0)
316 {
317 K = 1.0E+307;
318 }
319 else if (x < 9.0)
320 {
321 G4double y = x * x;
322 G4double C = -std::log(x/2.0) - 0.5772156649015329;
323 G4double q = 1.0;
324 G4double t = 0.0;
325 for (G4int i=1; i<51; i++)
326 {
327 q *= 0.25 * y / i / i;
328 t += 1.0 / i ;
329 K += q * (t+C);
330 }
331 K += C;
332 }
333 else
334 {
335 G4double y = 1.0 / x / x;
336 K = 0.5 / x / pI0(x) *
337 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7))))))));
338 }
339
340 return K;
341}
342////////////////////////////////////////////////////////////////////////////////
343//
344G4double G4Bessel::pK1 (G4double x)
345{
346 G4double K = 0.0;
347 if (x == 0.0)
348 K = 1.0E+307;
349 else
350 K = (1.0/x - pI1(x)*pK0(x)) / pI0(x);
351 return K;
352}
353////////////////////////////////////////////////////////////////////////////////
354//
A11
#define A11
Definition: G4ErrorMatrix.cc:825
A10
#define A10
Definition: G4ErrorMatrix.cc:824
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4Bessel::~G4Bessel
~G4Bessel()
Definition: G4Bessel.cc:69
G4Bessel::K0
G4double K0(G4double)
Definition: G4Bessel.cc:109
G4Bessel::I0
G4double I0(G4double)
Definition: G4Bessel.cc:73
G4Bessel::I1
G4double I1(G4double)
Definition: G4Bessel.cc:143
G4Bessel::K1
G4double K1(G4double)
Definition: G4Bessel.cc:181
G4Bessel::pI1
G4double pI1(G4double)
Definition: G4Bessel.cc:258
G4Bessel::pK0
G4double pK0(G4double)
Definition: G4Bessel.cc:303
G4Bessel::pI0
G4double pI0(G4double)
Definition: G4Bessel.cc:215
G4Bessel::pK1
G4double pK1(G4double)
Definition: G4Bessel.cc:344
G4Bessel::G4Bessel
G4Bessel()
Definition: G4Bessel.cc:65