Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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
42G4QBesIKJY::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
114G4QBesIKJY::~G4QBesIKJY(){;}
115
116G4double 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}
G4QBIType
G4QBIType
Definition: G4QBesIKJY.hh:43
BessJ0
@ BessJ0
Definition: G4QBesIKJY.hh:43
BessY0
@ BessY0
Definition: G4QBesIKJY.hh:44
BessJ1
@ BessJ1
Definition: G4QBesIKJY.hh:43
BessY1
@ BessY1
Definition: G4QBesIKJY.hh:44
BessK1
@ BessK1
Definition: G4QBesIKJY.hh:44
EBessK1
@ EBessK1
Definition: G4QBesIKJY.hh:44
EBessK0
@ EBessK0
Definition: G4QBesIKJY.hh:44
BessK0
@ BessK0
Definition: G4QBesIKJY.hh:44
EBessI1
@ EBessI1
Definition: G4QBesIKJY.hh:43
EBessI0
@ EBessI0
Definition: G4QBesIKJY.hh:43
BessI1
@ BessI1
Definition: G4QBesIKJY.hh:43
BessI0
@ BessI0
Definition: G4QBesIKJY.hh:43
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:69
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4QBesIKJY::operator()
G4double operator()(G4double x) const
Definition: G4QBesIKJY.cc:116
G4QBesIKJY::G4QBesIKJY
G4QBesIKJY(G4QBIType type=BessI0)
Definition: G4QBesIKJY.cc:42
V1
#define V1(a, b, c)
V2
#define V2(a, b, c)