Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/util/src/G4Bessel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 72 G4Bessel::G4Bessel ()
 73 {;}
 74 ////////////////////////////////////////////////////////////////////////////////
 75 //
 76 G4Bessel::~G4Bessel ()
 77 {;}
 78 ////////////////////////////////////////////////////////////////////////////////
 79 //
 80 G4double 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 //
116 G4double 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 //
150 G4double 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 //
188 G4double 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 //
222 G4double 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 //
265 G4double 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 //
310 G4double 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 //
351 G4double 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 //
362