Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QBesIKJY.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 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 //
29 // ---------------- G4QBesIKJY ----------------
30 // by Mikhail Kossov, Sept 1999.
31 // class for Bessel I0/I1 and K0/K1 functions in CHIPS Model
32 // -------------------------------------------------------------------
33 // Short description: Bessel functions class (can be substituted)
34 // -------------------------------------------------------------------
35 
36 //#define debug
37 //#define pdebug
38 
39 #include "G4QBesIKJY.hh"
40 
41 // Constructor
42 G4QBesIKJY::G4QBesIKJY(G4QBIType type)
43 {
44  ex=false;
45  switch (type)
46  {
47  case BessI0:
48  nu=0;
49  ik=true;
50  ij=true;
51  break;
52  case BessI1:
53  nu=1;
54  ik=true;
55  ij=true;
56  break;
57  case EBessI0:
58  nu=0;
59  ex=true;
60  ik=true;
61  ij=true;
62  break;
63  case EBessI1:
64  nu=1;
65  ex=true;
66  ik=true;
67  ij=true;
68  break;
69  case BessJ0:
70  nu=0;
71  ik=true;
72  ij=false;
73  break;
74  case BessJ1:
75  nu=1;
76  ik=true;
77  ij=false;
78  break;
79  case BessK0:
80  nu=0;
81  ik=false;
82  ij=true;
83  break;
84  case BessK1:
85  nu=1;
86  ik=false;
87  ij=true;
88  break;
89  case EBessK0:
90  nu=0;
91  ex=true;
92  ik=false;
93  ij=true;
94  break;
95  case EBessK1:
96  nu=1;
97  ex=true;
98  ik=false;
99  ij=true;
100  break;
101  case BessY0:
102  nu=0;
103  ik=false;
104  ij=false;
105  break;
106  case BessY1:
107  nu=1;
108  ik=false;
109  ij=false;
110  break;
111  }
112 }
113 
114 G4QBesIKJY::~G4QBesIKJY(){;}
115 
116 G4double G4QBesIKJY::operator() (G4double X) const
117 {
118  static const G4int nat1 = 15; // a # of attempts to reach the X<1 accuracy
119  static const G4int nat2 = nat1+nat1; // a # of attempts to reach the X<5 accuracy
120  static const G4int npi = 25;
121  static const G4int npil = npi-1;
122  static const G4int npk = 17;
123  static const G4int npkl = npk-1;
124  static const G4int npj = 20;
125  static const G4int npjl = npj-1;
126  static const G4complex CI(0,1);
127  static const G4double EPS = 1.E-15;
128  static const G4double Z1 = 1.;
129  static const G4double HF = Z1/2;
130  static const G4double R8 = HF/4;
131  static const G4double R32 = R8/4;
132  static const G4double PI = 3.14159265358979324;
133  static const G4double CE = 0.57721566490153286;
134  static const G4double PIH = PI/2;
135  static const G4double PI4 = PIH/2; // PI/4
136  static const G4double PI3 = PIH+PI4; // 3*PI/4
137  static const G4double RPIH = 2./PI;
138  static const G4double RPI2 = RPIH/4;
139 
140  static const G4double CI0[npi]={+1.00829205458740032E0,
141  +.00845122624920943E0,+.00012700630777567E0,+.00007247591099959E0,
142  +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0,
143  +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0,
144  -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0,
145  +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0,
146  -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0,
147  -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0,
148  +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0};
149 
150  static const G4double CI1[npi]={+.975800602326285926E0,
151  -.024467442963276385E0,-.000277205360763829E0,-.000009732146728020E0,
152  -.000000629724238640E0,-.000000065961142154E0,-.000000009613872919E0,
153  -.000000001401140901E0,-.000000000047563167E0,+.000000000081530681E0,
154  +.000000000035408148E0,+.000000000005102564E0,-.000000000001804409E0,
155  -.000000000001023594E0,-.000000000000052678E0,+.000000000000107094E0,
156  +.000000000000026120E0,-.000000000000009561E0,-.000000000000004713E0,
157  +.000000000000000829E0,+.000000000000000743E0,-.000000000000000080E0,
158  -.000000000000000117E0,+.000000000000000011E0,+.000000000000000019E0};
159 
160  static const G4double CK0[npk]={+.988408174230825800E0,-.011310504646928281E0,
161  +.000269532612762724E0,-.000011106685196665E0,+.000000632575108500E0,
162  -.000000045047337641E0,+.000000003792996456E0,-.000000000364547179E0,
163  +.000000000039043756E0,-.000000000004579936E0,+.000000000000580811E0,
164  -.000000000000078832E0,+.000000000000011360E0,-.000000000000001727E0,
165  +.000000000000000275E0,-.000000000000000046E0,+.000000000000000008E0};
166 
167  static const G4double CK1[npk]={+1.035950858772358331E0,+.035465291243331114E0,
168  -.000468475028166889E0,+.000016185063810053E0,-.000000845172048124E0,
169  +.000000057132218103E0,-.000000004645554607E0,+.000000000435417339E0,
170  -.000000000045757297E0,+.000000000005288133E0,-.000000000000662613E0,
171  +.000000000000089048E0,-.000000000000012726E0,+.000000000000001921E0,
172  -.000000000000000305E0,+.000000000000000050E0,-.000000000000000009E0};
173 
174  static const G4double CA[npk]={+.157727971474890120E0,-.008723442352852221E0,
175  +.265178613203336810E0,-.370094993872649779E0,+.158067102332097261E0,
176  -.034893769411408885E0,+.004819180069467605E0,-.000460626166206275E0,
177  +.000032460328821005E0,-.000001761946907762E0,+.000000076081635924E0,
178  -.000000002679253531E0,+.000000000078486963E0,-.000000000001943835E0,
179  +.000000000000041253E0,-.000000000000000759E0,+.000000000000000012E0};
180 
181  static const G4double CB[npk]={-.021505111449657551E0,-.275118133043518791E0,
182  +.198605634702554156E0,+.234252746109021802E0,-.165635981713650413E0,
183  +.044621379540669282E0,-.006932286291523188E0,+.000719117403752303E0,
184  -.000053925079722939E0,+.000003076493288108E0,-.000000138457181230E0,
185  +.000000005051054369E0,-.000000000152582850E0,+.000000000003882867E0,
186  -.000000000000084429E0,+.000000000000001587E0,-.000000000000000026E0};
187 
188  static const G4complex CC[npj]={
189  G4complex(+.998988089858965153E0,-.012331520578544144E0),
190  G4complex(-.001338428549971856E0,-.012249496281259475E0),
191  G4complex(-.000318789878061893E0,+.000096494184993423E0),
192  G4complex(+.000008511232210657E0,+.000013655570490357E0),
193  G4complex(+.000000691542349139E0,-.000000851806644426E0),
194  G4complex(-.000000090770101537E0,-.000000027244053414E0),
195  G4complex(+.000000001454928079E0,+.000000009646421338E0),
196  G4complex(+.000000000926762487E0,-.000000000683347518E0),
197  G4complex(-.000000000139166198E0,-.000000000060627380E0),
198  G4complex(+.000000000003237975E0,+.000000000021695716E0),
199  G4complex(+.000000000002535357E0,-.000000000002304899E0),
200  G4complex(-.000000000000559090E0,-.000000000000122554E0),
201  G4complex(+.000000000000041919E0,+.000000000000092314E0),
202  G4complex(+.000000000000008733E0,-.000000000000016778E0),
203  G4complex(-.000000000000003619E0,+.000000000000000754E0),
204  G4complex(+.000000000000000594E0,+.000000000000000462E0),
205  G4complex(-.000000000000000010E0,-.000000000000000159E0),
206  G4complex(-.000000000000000024E0,+.000000000000000025E0),
207  G4complex(+.000000000000000008E0,+.000000000000000000E0),
208  G4complex(-.000000000000000001E0,-.000000000000000001E0)};
209 
210  static const G4double CD[npk]={+0.648358770605264921E0,-1.191801160541216873E0,
211  +1.287994098857677620E0,-0.661443934134543253E0,+0.177709117239728283E0,
212  -0.029175524806154208E0,+0.003240270182683857E0,-0.000260444389348581E0,
213  +0.000015887019239932E0,-0.000000761758780540E0,+0.000000029497070073E0,
214  -0.000000000942421298E0,+0.000000000025281237E0,-0.000000000000577740E0,
215  +0.000000000000011386E0,-0.000000000000000196E0,+0.000000000000000003E0};
216 
217  static const G4double EE[npk]={-.040172946544414076E0,-.444447147630558063E0,
218  -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0,
219  +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0,
220  -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0,
221  +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0,
222  -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0};
223 
224  static const G4complex CF[npj]={
225  G4complex(+1.001702234853820996E0,+.037261715000537654E0),
226  G4complex(+.002255572846561180E0,+.037145322479807690E0),
227  G4complex(+.000543216487508013E0,-.000137263238201907E0),
228  G4complex(-.000011179461895408E0,-.000019851294687597E0),
229  G4complex(-.000000946901382392E0,+.000001070014057386E0),
230  G4complex(+.000000111032677121E0,+.000000038305261714E0),
231  G4complex(-.000000001294398927E0,-.000000011628723277E0),
232  G4complex(-.000000001114905944E0,+.000000000759733092E0),
233  G4complex(+.000000000157637232E0,+.000000000075476075E0),
234  G4complex(-.000000000002830457E0,-.000000000024752781E0),
235  G4complex(-.000000000002932169E0,+.000000000002493893E0),
236  G4complex(+.000000000000617809E0,+.000000000000156198E0),
237  G4complex(-.000000000000043162E0,-.000000000000103385E0),
238  G4complex(-.000000000000010133E0,+.000000000000018129E0),
239  G4complex(+.000000000000003973E0,-.000000000000000709E0),
240  G4complex(-.000000000000000632E0,-.000000000000000520E0),
241  G4complex(+.000000000000000006E0,+.000000000000000172E0),
242  G4complex(+.000000000000000027E0,-.000000000000000026E0),
243  G4complex(-.000000000000000008E0,-.000000000000000000E0),
244  G4complex(+.000000000000000001E0,+.000000000000000001E0)};
245  // -------------------------------------------------------------------------------------
246  G4double H=0.; // Prototype of the result
247  if (ij) // I/K Bessel functions
248  {
249  if (ik) // I0/I1/EI0/EI1 Bessel functions (symmetric)
250  {
251  G4double V=std::abs(X);
252  G4double CJ=0.; // Prototype of the element of the CI0/CI1 matrix
253  if (V < 8.)
254  {
255  G4double HFV=HF*V;
256  G4double Y=HFV*HFV;
257  G4int V3=nu+1;
258  G4int XL=V3+1;
259  G4int XLI=XL+1;
260  G4int XLD=XLI+1;
261  G4int W1=XLD+1;
262  G4double A0=1.;
263  G4double DY=Y+Y;
264  G4double A1=1.+DY/(XLI*V3);
265  G4double A2=1.+Y*(4.+(DY+Y)/(XLD*XL))/(W1*V3);
266  G4double B0=1.;
267  G4double B1=1.-Y/XLI;
268  G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1;
269  G4int V1=3-XL;
270  G4double V2=V3+V3;
271  G4double C=0.;
272  for (G4int N=3; N<=30; N++)
273  {
274  G4double C0=C;
275  G4double FN=N;
276  W1=W1+2;
277  G4int W2=W1-1;
278  G4int W3=W2-1;
279  G4int W4=W3-1;
280  G4int W5=W4-1;
281  G4int W6=W5-1;
282  V1=V1+1;
283  V2=V2+1;
284  V3=V3+1;
285  G4double U1=FN*W4;
286  G4double E=V3/(U1*W3);
287  G4double U2=E*Y;
288  G4double F1=1.+Y*V1/(U1*W1);
289  G4double F2=(1.+Y*V2/(V3*W2*W5))*U2;
290  G4double F3=-Y*Y*U2/(W4*W5*W5*W6);
291  G4double A=F1*A2+F2*A1+F3*A0;
292  G4double B=F1*B2+F2*B1+F3*B0;
293  C=A/B;
294  if (std::abs(C0-C) < EPS*std::abs(C)) break;
295  A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
296  }
297  H=C;
298  if (nu==1) H*=HF*X;
299  if (ex) H*=std::exp(-V);
300  }
301  else
302  {
303  G4double P=16./V-1.;
304  G4double ALFA=P+P;
305  G4double B1=0.;
306  G4double B2=0.;
307  for (G4int I=npil; I>=0; I--)
308  {
309  if (!nu) CJ=CI0[I];
310  else CJ=CI1[I];
311  G4double B0=CJ+ALFA*B1-B2;
312  B2=B1;
313  B1=B0;
314  }
315  H=std::sqrt(RPI2/V)*(B1-P*B2);
316  if (nu && X < 0.) H=-H;
317  if (!ex) H*=std::exp(V);
318  }
319  }
320  else // K0/K1/EK0/EK1 Bessel functions
321  {
322 #ifdef debug
323  G4cout<<"G4BesIKJY: ------------>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
324 #endif
325  G4double CK=0.; // Prototype of the element of the CI0/CI1 matrix
326  if (X < 0.)
327  {
328  G4cout<<"G4BesIKJY::NegativeArg in K-BessFun X="<<X<<", n="<<nu<<",E="<<ex<<G4endl;
329  return H;
330  }
331  else if (X < 1.)
332  {
333 #ifdef debug
334  G4cout<<"G4BesIKJY: -->> [ X < 1 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
335 #endif
336  G4double B=HF*X;
337  G4double BK=-(std::log(B)+CE);
338  G4double F=BK;
339  G4double P=HF;
340  G4double Q=HF;
341  G4double C=1.;
342  G4double D=B*B;
343  G4double BK1=P;
344  for (G4int N=1; N<=nat1; N++) // @@ "nat" can be increased
345  {
346  G4double FN=N;
347  P/=FN;
348  Q/=FN;
349  F=(F+P+Q)/FN;
350  C*=D/FN;
351  G4double G=C*(P-FN*F);
352  G4double R=C*F;
353  BK=BK+R;
354  BK1=BK1+G;
355  if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break;
356  }
357  if (nu==1) H=BK1/B;
358  else H=BK;
359  if (ex) H*=std::exp(X);
360  }
361  else if (X < 5.)
362  {
363 #ifdef debug
364  G4cout<<"G4BesIKJY: -->> [ X < 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
365 #endif
366  G4int NUS=0; // @@ NU**2 for future NU>1 applications
367  if (nu==1) NUS=1;
368  G4double DNUS=NUS+NUS;
369  G4double XN=DNUS+DNUS;
370  G4double A=9.-XN;
371  G4double B=25.-XN;
372  G4double C=768*X*X;
373  G4double HX=16*X;
374  G4double C0=HX+HX+HX;;
375  G4double A0=1.;
376  G4double A1=(HX+7.+XN)/A;
377  G4double A2=(C+C0*(XN+23.)+XN*(XN+62.)+129.)/(A*B);
378  G4double B0=1.;
379  G4double B1=(HX+9.-XN)/A;
380  G4double B2=(C+C0*B)/(A*B)+1.;
381  C=0.;
382  for (G4int N=3; N<=nat2; N++)
383  {
384  C0=C;
385  G4double FN=N;
386  G4double FN2=FN+FN;
387  G4double FNP=FN2+1.;
388  G4double FN1=FN2-1.;
389  G4double FNM=FN1-4.;
390  G4double FN3=FN1/(FN2-3.);
391  G4double FN4=12*FN*FN-(1.-XN);
392  G4double FN5=16*FN1*X;
393  G4double RAN=1./(FNP*FNP-XN);
394  G4double F1=FN3*(FN4-20*FN)+FN5;
395  G4double F2=28*FN-FN4-8.+FN5;
396  G4double F3=FN3*(FNM*FNM-XN);
397  A=(F1*A2+F2*A1+F3*A0)*RAN;
398  B=(F1*B2+F2*B1+F3*B0)*RAN;
399  C=A/B;
400  if (std::abs(C0-C) < EPS*std::abs(C)) break;
401  A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B;
402  }
403  H=C/std::sqrt(RPIH*X);
404  if (!ex) H*=std::exp(-X);
405  }
406  else
407  {
408 #ifdef debug
409  G4cout<<"G4BesIKJY: ->> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
410 #endif
411  G4double P=10./X-1.;
412  G4double ALFA=P+P;
413  G4double B1=0.;
414  G4double B2=0.;
415 #ifdef debug
416  G4cout<<"G4BesIKJY: ->> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
417 #endif
418  for (G4int I=npkl; I>=0; I--)
419  {
420  if (!nu) CK=CK0[I];
421  else CK=CK1[I];
422  G4double B0=CK+ALFA*B1-B2;
423  B2=B1;
424  B1=B0;
425  }
426  H=std::sqrt(PIH/X)*(B1-P*B2);
427  if (!ex) H*=std::exp(-X);
428  }
429  }
430  }
431  else
432  {
433  if (!ik && X < 0.)
434  {
435  G4cout<<"G4BesIKJY::NegativeArgument in Y BesselFunction X="<<X<<", nu="<<nu<<G4endl;
436  return H;
437  }
438  G4double V=std::abs(X);
439  if (!nu) // J0/Y0 Bessel functions
440  {
441  if (V < 8.)
442  {
443  G4double P=R32*V*V-1.;
444  G4double ALFA=P+P;
445  G4double B1=0.;
446  G4double B2=0.;
447  for (G4int IT=npkl; IT>=0; IT--)
448  {
449  G4double B0=CA[IT]+ALFA*B1-B2;
450  B2=B1;
451  B1=B0;
452  }
453  H=B1-P*B2;
454  if (!ik)
455  {
456  B1=0.;
457  B2=0.;
458  for (G4int JT=npkl; JT>=0; JT--)
459  {
460  G4double B0=CB[JT]+ALFA*B1-B2;
461  B2=B1;
462  B1=B0;
463  }
464  H=RPIH*(CE+std::log(HF*X))*H+B1-P*B2;
465  }
466  }
467  else
468  {
469  G4double P=10./V-1.;
470  G4double ALFA=P+P;
471  G4complex CB1(0.,0.);
472  G4complex CB2(0.,0.);
473  for (G4int IT=npjl; IT>=0; IT--)
474  {
475  G4complex CB0=CC[IT]+ALFA*CB1-CB2;
476  CB2=CB1;
477  CB1=CB0;
478  }
479  CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2);
480  if (ik) H=real(CB1);
481  else H=real(-CI*CB1);
482  }
483  }
484  else // J1/Y1 Bessel functions
485  {
486  if (V < 8.)
487  {
488  G4double Y=R8*V;
489  G4double Y2=Y*Y;
490  G4double P=Y2+Y2-1.;
491  G4double ALFA=P+P;
492  G4double B1=0.;
493  G4double B2=0.;
494  for (G4int IT=npkl; IT>=0; IT--)
495  {
496  G4double B0=CD[IT]+ALFA*B1-B2;
497  B2=B1;
498  B1=B0;
499  }
500  H=Y*(B1-P*B2);
501  if (!ik)
502  {
503  B1=0.;
504  B2=0.;
505  for (G4int JT=npkl; JT>=0; JT--)
506  {
507  G4double B0=EE[JT]+ALFA*B1-B2;
508  B2=B1;
509  B1=B0;
510  }
511  H=RPIH*((CE+std::log(HF*X))*H-1./X)+Y*(B1-B2);
512  }
513  }
514  else
515  {
516  G4double P=10./V-1.;
517  G4double ALFA=P+P;
518  G4complex CB1(0.,0.);
519  G4complex CB2(0.,0.);
520  for (G4int IT=npjl; IT>=0; IT--)
521  {
522  G4complex CB0=CF[IT]+ALFA*CB1-CB2;
523  CB2=CB1;
524  CB1=CB0;
525  }
526  CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI3))*(CB1-P*CB2);
527  if (ik) H=real(CB1);
528  else H=real(-CI*CB1);
529  }
530  if (X < 0.) H=-H;
531  }
532  }
533  return H;
534 }
Generated on Sat May 25 2013 14:33:45 for Geant4 by   doxygen 1.8.4