source: trunk/source/processes/hadronic/util/src/G4Bessel.cc@ 1055

Last change on this file since 1055 was 819, checked in by garnier, 17 years ago

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