Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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// 06 August 2015, A. Ribon, CERN
59// Migrated to G4Exp, G4Log and G4Pow.
60//
61// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62////////////////////////////////////////////////////////////////////////////////
63//
64#include "G4Bessel.hh"
65#include "G4PhysicalConstants.hh"
66
67#include "G4Exp.hh"
68#include "G4Log.hh"
69#include "G4Pow.hh"
70////////////////////////////////////////////////////////////////////////////////
71//
72G4Bessel::G4Bessel ()
73{;}
74////////////////////////////////////////////////////////////////////////////////
75//
76G4Bessel::~G4Bessel ()
77{;}
78////////////////////////////////////////////////////////////////////////////////
79//
80G4double G4Bessel::I0 (G4double x)
81{
82 const G4double P1 = 1.0,
83 P2 = 3.5156229,
84 P3 = 3.0899424,
85 P4 = 1.2067492,
86 P5 = 0.2659732,
87 P6 = 0.0360768,
88 P7 = 0.0045813;
89 const G4double Q1 = 0.39894228,
90 Q2 = 0.01328592,
91 Q3 = 0.00225319,
92 Q4 =-0.00157565,
93 Q5 = 0.00916281,
94 Q6 =-0.02057706,
95 Q7 = 0.02635537,
96 Q8 =-0.01647633,
97 Q9 = 0.00392377;
98
99 G4double I = 0.0;
100 if (std::abs(x) < 3.75)
101 {
102 G4double y = G4Pow::GetInstance()->powN(x/3.75, 2);
103 I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
104 }
105 else
106 {
107 G4double ax = std::abs(x);
108 G4double y = 3.75/ax;
109 I = G4Exp(ax) / std::sqrt(ax) *
110 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
111 }
112 return I;
113}
114////////////////////////////////////////////////////////////////////////////////
115//
116G4double G4Bessel::K0 (G4double x)
117{
118 const G4double P1 =-0.57721566,
119 P2 = 0.42278420,
120 P3 = 0.23069756,
121 P4 = 0.03488590,
122 P5 = 0.00262698,
123 P6 = 0.00010750,
124 P7 = 0.00000740;
125 const G4double Q1 = 1.25331414,
126 Q2 =-0.07832358,
127 Q3 = 0.02189568,
128 Q4 =-0.01062446,
129 Q5 = 0.00587872,
130 Q6 =-0.00251540,
131 Q7 = 0.00053208;
132
133 G4double K = 0.0;
134 if (x <= 2.0)
135 {
136 G4double y = x * x / 4.0;
137 K = (-G4Log(x/2.0)) * I0(x) +
138 P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))));
139 }
140 else
141 {
142 G4double y = 2.0 / x;
143 K = G4Exp(-x) / std::sqrt(x) *
144 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
145 }
146 return K;
147}
148////////////////////////////////////////////////////////////////////////////////
149//
150G4double G4Bessel::I1 (G4double x)
151{
152 const G4double P1 = 0.5,
153 P2 = 0.87890594,
154 P3 = 0.51498869,
155 P4 = 0.15084934,
156 P5 = 0.02658733,
157 P6 = 0.00301532,
158 P7 = 0.00032411;
159 const G4double Q1 = 0.39894228,
160 Q2 =-0.03988024,
161 Q3 =-0.00362018,
162 Q4 = 0.00163801,
163 Q5 =-0.01031555,
164 Q6 = 0.02282967,
165 Q7 =-0.02895312,
166 Q8 = 0.01787654,
167 Q9 =-0.00420059;
168
169 G4double I = 0.0;
170 if (std::abs(x) < 3.75)
171 {
172 G4double ax = std::abs(x);
173 G4double y = G4Pow::GetInstance()->powN(x/3.75, 2);
174 I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
175 }
176 else
177 {
178 G4double ax = std::abs(x);
179 G4double y = 3.75/ax;
180 I = G4Exp(ax) / std::sqrt(ax) *
181 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9))))))));
182 }
183 if (x < 0.0) I = -I;
184 return I;
185}
186////////////////////////////////////////////////////////////////////////////////
187//
188G4double G4Bessel::K1 (G4double x)
189{
190 const G4double P1 = 1.0,
191 P2 = 0.15443144,
192 P3 =-0.67278579,
193 P4 =-0.18156897,
194 P5 =-0.01919402,
195 P6 =-0.00110404,
196 P7 =-0.00004686;
197 const G4double Q1 = 1.25331414,
198 Q2 = 0.23498619,
199 Q3 =-0.03655620,
200 Q4 = 0.01504268,
201 Q5 =-0.00780353,
202 Q6 = 0.00325614,
203 Q7 =-0.00068245;
204
205 G4double K = 0.0;
206 if (x <= 2.0)
207 {
208 G4double y = x * x / 4.0;
209 K = G4Log(x/2.0)*I1(x) + 1.0/x *
210 (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))));
211 }
212 else
213 {
214 G4double y = 2.0 / x;
215 K = G4Exp(-x) / std::sqrt(x) *
216 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7))))));
217 }
218 return K;
219}
220////////////////////////////////////////////////////////////////////////////////
221//
222G4double G4Bessel::pI0 (G4double x)
223{
224 const G4double A0 = 0.1250000000000E+00,
225 A1 = 7.0312500000000E-02,
226 A2 = 7.3242187500000E-02,
227 A3 = 1.1215209960938E-01,
228 A4 = 2.2710800170898E-01,
229 A5 = 5.7250142097473E-01,
230 A6 = 1.7277275025845E+00,
231 A7 = 6.0740420012735E+00,
232 A8 = 2.4380529699556E+01,
233 A9 = 1.1001714026925E+02,
234 A10 = 5.5133589612202E+02,
235 A11 = 3.0380905109224E+03;
236
237 G4double I = 0.0;
238 if (x == 0.0)
239 {
240 I = 1.0;
241 }
242 else if (x < 18.0)
243 {
244 I = 1.0;
245 G4double y = x * x;
246 G4double q = 1.0;
247 for (G4int i=1; i<101; i++)
248 {
249 q *= 0.25 * y / i / i;
250 I += q;
251 if (std::abs(q/I) < 1.0E-15) break;
252 }
253 }
254 else
255 {
256 G4double y = 1.0 / x;
257 I = G4Exp(x) / std::sqrt(twopi*x) *
258 (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))))))))))));
259 }
260
261 return I;
262}
263////////////////////////////////////////////////////////////////////////////////
264//
265G4double G4Bessel::pI1 (G4double x)
266{
267 const G4double A0 = -0.3750000000000E+00,
268 A1 = -1.1718750000000E-01,
269 A2 = -1.0253906250000E-01,
270 A3 = -1.4419555664063E-01,
271 A4 = -2.775764465332E-01,
272 A5 = -6.7659258842468E-01,
273 A6 = -1.9935317337513E+00,
274 A7 = -6.8839142681099E+00,
275 A8 = -2.7248827311269E+01,
276 A9 = -1.2159789187654E+02,
277 A10 = -6.0384407670507E+02,
278 A11 = -3.3022722944809E+03;
279
280 G4double I = 0.0;
281 if (x == 0.0)
282 {
283 I = 0.0;
284 }
285 else if (x < 18.0)
286 {
287 I = 1.0;
288 G4double y = x * x;
289 G4double q = 1.0;
290 for (G4int i=1; i<101; i++)
291 {
292 q *= 0.25 * y / i / (i+1.0);
293 I += q;
294 if (std::abs(q/I) < 1.0E-15) break;
295 }
296 I *= 0.5 * x;
297
298 }
299 else
300 {
301 G4double y = 1.0 / x;
302 I = G4Exp(x) / std::sqrt(twopi*x) *
303 (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))))))))))));
304 }
305
306 return I;
307}
308////////////////////////////////////////////////////////////////////////////////
309//
310G4double G4Bessel::pK0 (G4double x)
311{
312 const G4double A0 = 0.1250000000000E+00,
313 A1 = 0.2109375000000E+00,
314 A2 = 1.0986328125000E+00,
315 A3 = 1.1775970458984E+01,
316 A4 = 2.1461706161499E+02,
317 A5 = 5.9511522710323E+03,
318 A6 = 2.3347645606175E+05,
319 A7 = 1.2312234987631E+07;
320
321 G4double K = 0.0;
322 if (x == 0.0)
323 {
324 K = 1.0E+307;
325 }
326 else if (x < 9.0)
327 {
328 G4double y = x * x;
329 G4double C = -G4Log(x/2.0) - 0.5772156649015329;
330 G4double q = 1.0;
331 G4double t = 0.0;
332 for (G4int i=1; i<51; i++)
333 {
334 q *= 0.25 * y / i / i;
335 t += 1.0 / i ;
336 K += q * (t+C);
337 }
338 K += C;
339 }
340 else
341 {
342 G4double y = 1.0 / x / x;
343 K = 0.5 / x / pI0(x) *
344 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7))))))));
345 }
346
347 return K;
348}
349////////////////////////////////////////////////////////////////////////////////
350//
351G4double G4Bessel::pK1 (G4double x)
352{
353 G4double K = 0.0;
354 if (x == 0.0)
355 K = 1.0E+307;
356 else
357 K = (1.0/x - pI1(x)*pK0(x)) / pI0(x);
358 return K;
359}
360////////////////////////////////////////////////////////////////////////////////
361//
A11
#define A11
Definition: G4ErrorMatrix.cc:875
A10
#define A10
Definition: G4ErrorMatrix.cc:874
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4Bessel::~G4Bessel
~G4Bessel()
Definition: G4Bessel.cc:76
G4Bessel::K0
G4double K0(G4double)
Definition: G4Bessel.cc:116
G4Bessel::I0
G4double I0(G4double)
Definition: G4Bessel.cc:80
G4Bessel::I1
G4double I1(G4double)
Definition: G4Bessel.cc:150
G4Bessel::K1
G4double K1(G4double)
Definition: G4Bessel.cc:188
G4Bessel::pI1
G4double pI1(G4double)
Definition: G4Bessel.cc:265
G4Bessel::pK0
G4double pK0(G4double)
Definition: G4Bessel.cc:310
G4Bessel::pI0
G4double pI0(G4double)
Definition: G4Bessel.cc:222
G4Bessel::pK1
G4double pK1(G4double)
Definition: G4Bessel.cc:351
G4Bessel::G4Bessel
G4Bessel()
Definition: G4Bessel.cc:72
G4Pow::GetInstance
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4Pow::powN
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:162