Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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