Geant4  10.02.p01
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 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 //
64 #include "G4Bessel.hh"
65 #include "G4PhysicalConstants.hh"
66 
67 #include "G4Exp.hh"
68 #include "G4Log.hh"
69 #include "G4Pow.hh"
71 //
73 {;}
75 //
77 {;}
79 //
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::fabs(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::fabs(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 }
115 //
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 }
149 //
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::fabs(x) < 3.75)
171  {
172  G4double ax = std::fabs(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::fabs(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 }
187 //
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 }
221 //
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 }
264 //
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 }
309 //
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 }
350 //
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 }
361 //
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double K1(G4double)
Definition: G4Bessel.cc:188
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
static const G4double * P1[nN]
G4double pI0(G4double)
Definition: G4Bessel.cc:222
~G4Bessel()
Definition: G4Bessel.cc:76
#define A10
G4double I0(G4double)
Definition: G4Bessel.cc:80
G4double pK1(G4double)
Definition: G4Bessel.cc:351
G4double pI1(G4double)
Definition: G4Bessel.cc:265
double C(double temp)
int G4int
Definition: G4Types.hh:78
G4double I1(G4double)
Definition: G4Bessel.cc:150
static const double twopi
Definition: G4SIunits.hh:75
G4double K0(G4double)
Definition: G4Bessel.cc:116
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double pK0(G4double)
Definition: G4Bessel.cc:310
const G4double x[NPOINTSGL]
#define A11
static const G4double * P2[nN]
double G4double
Definition: G4Types.hh:76
G4Bessel()
Definition: G4Bessel.cc:72