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