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

Last change on this file since 1183 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 10.5 KB
Line 
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.