Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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)