#include <G4QBesIKJY.hh>
Public Member Functions | |
| G4QBesIKJY (G4QBIType type=BessI0) | |
| ~G4QBesIKJY () | |
| G4double | operator() (G4double x) const |
Definition at line 45 of file G4QBesIKJY.hh.
| G4QBesIKJY::G4QBesIKJY | ( | G4QBIType | type = BessI0 |
) |
Definition at line 42 of file G4QBesIKJY.cc.
References BessI0, BessI1, BessJ0, BessJ1, BessK0, BessK1, BessY0, BessY1, EBessI0, EBessI1, EBessK0, and EBessK1.
00043 { 00044 ex=false; 00045 switch (type) 00046 { 00047 case BessI0: 00048 nu=0; 00049 ik=true; 00050 ij=true; 00051 break; 00052 case BessI1: 00053 nu=1; 00054 ik=true; 00055 ij=true; 00056 break; 00057 case EBessI0: 00058 nu=0; 00059 ex=true; 00060 ik=true; 00061 ij=true; 00062 break; 00063 case EBessI1: 00064 nu=1; 00065 ex=true; 00066 ik=true; 00067 ij=true; 00068 break; 00069 case BessJ0: 00070 nu=0; 00071 ik=true; 00072 ij=false; 00073 break; 00074 case BessJ1: 00075 nu=1; 00076 ik=true; 00077 ij=false; 00078 break; 00079 case BessK0: 00080 nu=0; 00081 ik=false; 00082 ij=true; 00083 break; 00084 case BessK1: 00085 nu=1; 00086 ik=false; 00087 ij=true; 00088 break; 00089 case EBessK0: 00090 nu=0; 00091 ex=true; 00092 ik=false; 00093 ij=true; 00094 break; 00095 case EBessK1: 00096 nu=1; 00097 ex=true; 00098 ik=false; 00099 ij=true; 00100 break; 00101 case BessY0: 00102 nu=0; 00103 ik=false; 00104 ij=false; 00105 break; 00106 case BessY1: 00107 nu=1; 00108 ik=false; 00109 ij=false; 00110 break; 00111 } 00112 }
| G4QBesIKJY::~G4QBesIKJY | ( | ) |
Definition at line 116 of file G4QBesIKJY.cc.
References G4cout, G4endl, IT, V1, and V2.
00117 { 00118 static const G4int nat1 = 15; // a # of attempts to reach the X<1 accuracy 00119 static const G4int nat2 = nat1+nat1; // a # of attempts to reach the X<5 accuracy 00120 static const G4int npi = 25; 00121 static const G4int npil = npi-1; 00122 static const G4int npk = 17; 00123 static const G4int npkl = npk-1; 00124 static const G4int npj = 20; 00125 static const G4int npjl = npj-1; 00126 static const G4complex CI(0,1); 00127 static const G4double EPS = 1.E-15; 00128 static const G4double Z1 = 1.; 00129 static const G4double HF = Z1/2; 00130 static const G4double R8 = HF/4; 00131 static const G4double R32 = R8/4; 00132 static const G4double PI = 3.14159265358979324; 00133 static const G4double CE = 0.57721566490153286; 00134 static const G4double PIH = PI/2; 00135 static const G4double PI4 = PIH/2; // PI/4 00136 static const G4double PI3 = PIH+PI4; // 3*PI/4 00137 static const G4double RPIH = 2./PI; 00138 static const G4double RPI2 = RPIH/4; 00139 00140 static const G4double CI0[npi]={+1.00829205458740032E0, 00141 +.00845122624920943E0,+.00012700630777567E0,+.00007247591099959E0, 00142 +.00000513587726878E0,+.00000056816965808E0,+.00000008513091223E0, 00143 +.00000001238425364E0,+.00000000029801672E0,-.00000000078956698E0, 00144 -.00000000033127128E0,-.00000000004497339E0,+.00000000001799790E0, 00145 +.00000000000965748E0,+.00000000000038604E0,-.00000000000104039E0, 00146 -.00000000000023950E0,+.00000000000009554E0,+.00000000000004443E0, 00147 -.00000000000000859E0,-.00000000000000709E0,+.00000000000000087E0, 00148 +.00000000000000112E0,-.00000000000000012E0,-.00000000000000018E0}; 00149 00150 static const G4double CI1[npi]={+.975800602326285926E0, 00151 -.024467442963276385E0,-.000277205360763829E0,-.000009732146728020E0, 00152 -.000000629724238640E0,-.000000065961142154E0,-.000000009613872919E0, 00153 -.000000001401140901E0,-.000000000047563167E0,+.000000000081530681E0, 00154 +.000000000035408148E0,+.000000000005102564E0,-.000000000001804409E0, 00155 -.000000000001023594E0,-.000000000000052678E0,+.000000000000107094E0, 00156 +.000000000000026120E0,-.000000000000009561E0,-.000000000000004713E0, 00157 +.000000000000000829E0,+.000000000000000743E0,-.000000000000000080E0, 00158 -.000000000000000117E0,+.000000000000000011E0,+.000000000000000019E0}; 00159 00160 static const G4double CK0[npk]={+.988408174230825800E0,-.011310504646928281E0, 00161 +.000269532612762724E0,-.000011106685196665E0,+.000000632575108500E0, 00162 -.000000045047337641E0,+.000000003792996456E0,-.000000000364547179E0, 00163 +.000000000039043756E0,-.000000000004579936E0,+.000000000000580811E0, 00164 -.000000000000078832E0,+.000000000000011360E0,-.000000000000001727E0, 00165 +.000000000000000275E0,-.000000000000000046E0,+.000000000000000008E0}; 00166 00167 static const G4double CK1[npk]={+1.035950858772358331E0,+.035465291243331114E0, 00168 -.000468475028166889E0,+.000016185063810053E0,-.000000845172048124E0, 00169 +.000000057132218103E0,-.000000004645554607E0,+.000000000435417339E0, 00170 -.000000000045757297E0,+.000000000005288133E0,-.000000000000662613E0, 00171 +.000000000000089048E0,-.000000000000012726E0,+.000000000000001921E0, 00172 -.000000000000000305E0,+.000000000000000050E0,-.000000000000000009E0}; 00173 00174 static const G4double CA[npk]={+.157727971474890120E0,-.008723442352852221E0, 00175 +.265178613203336810E0,-.370094993872649779E0,+.158067102332097261E0, 00176 -.034893769411408885E0,+.004819180069467605E0,-.000460626166206275E0, 00177 +.000032460328821005E0,-.000001761946907762E0,+.000000076081635924E0, 00178 -.000000002679253531E0,+.000000000078486963E0,-.000000000001943835E0, 00179 +.000000000000041253E0,-.000000000000000759E0,+.000000000000000012E0}; 00180 00181 static const G4double CB[npk]={-.021505111449657551E0,-.275118133043518791E0, 00182 +.198605634702554156E0,+.234252746109021802E0,-.165635981713650413E0, 00183 +.044621379540669282E0,-.006932286291523188E0,+.000719117403752303E0, 00184 -.000053925079722939E0,+.000003076493288108E0,-.000000138457181230E0, 00185 +.000000005051054369E0,-.000000000152582850E0,+.000000000003882867E0, 00186 -.000000000000084429E0,+.000000000000001587E0,-.000000000000000026E0}; 00187 00188 static const G4complex CC[npj]={ 00189 G4complex(+.998988089858965153E0,-.012331520578544144E0), 00190 G4complex(-.001338428549971856E0,-.012249496281259475E0), 00191 G4complex(-.000318789878061893E0,+.000096494184993423E0), 00192 G4complex(+.000008511232210657E0,+.000013655570490357E0), 00193 G4complex(+.000000691542349139E0,-.000000851806644426E0), 00194 G4complex(-.000000090770101537E0,-.000000027244053414E0), 00195 G4complex(+.000000001454928079E0,+.000000009646421338E0), 00196 G4complex(+.000000000926762487E0,-.000000000683347518E0), 00197 G4complex(-.000000000139166198E0,-.000000000060627380E0), 00198 G4complex(+.000000000003237975E0,+.000000000021695716E0), 00199 G4complex(+.000000000002535357E0,-.000000000002304899E0), 00200 G4complex(-.000000000000559090E0,-.000000000000122554E0), 00201 G4complex(+.000000000000041919E0,+.000000000000092314E0), 00202 G4complex(+.000000000000008733E0,-.000000000000016778E0), 00203 G4complex(-.000000000000003619E0,+.000000000000000754E0), 00204 G4complex(+.000000000000000594E0,+.000000000000000462E0), 00205 G4complex(-.000000000000000010E0,-.000000000000000159E0), 00206 G4complex(-.000000000000000024E0,+.000000000000000025E0), 00207 G4complex(+.000000000000000008E0,+.000000000000000000E0), 00208 G4complex(-.000000000000000001E0,-.000000000000000001E0)}; 00209 00210 static const G4double CD[npk]={+0.648358770605264921E0,-1.191801160541216873E0, 00211 +1.287994098857677620E0,-0.661443934134543253E0,+0.177709117239728283E0, 00212 -0.029175524806154208E0,+0.003240270182683857E0,-0.000260444389348581E0, 00213 +0.000015887019239932E0,-0.000000761758780540E0,+0.000000029497070073E0, 00214 -0.000000000942421298E0,+0.000000000025281237E0,-0.000000000000577740E0, 00215 +0.000000000000011386E0,-0.000000000000000196E0,+0.000000000000000003E0}; 00216 00217 static const G4double EE[npk]={-.040172946544414076E0,-.444447147630558063E0, 00218 -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0, 00219 +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0, 00220 -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0, 00221 +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0, 00222 -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0}; 00223 00224 static const G4complex CF[npj]={ 00225 G4complex(+1.001702234853820996E0,+.037261715000537654E0), 00226 G4complex(+.002255572846561180E0,+.037145322479807690E0), 00227 G4complex(+.000543216487508013E0,-.000137263238201907E0), 00228 G4complex(-.000011179461895408E0,-.000019851294687597E0), 00229 G4complex(-.000000946901382392E0,+.000001070014057386E0), 00230 G4complex(+.000000111032677121E0,+.000000038305261714E0), 00231 G4complex(-.000000001294398927E0,-.000000011628723277E0), 00232 G4complex(-.000000001114905944E0,+.000000000759733092E0), 00233 G4complex(+.000000000157637232E0,+.000000000075476075E0), 00234 G4complex(-.000000000002830457E0,-.000000000024752781E0), 00235 G4complex(-.000000000002932169E0,+.000000000002493893E0), 00236 G4complex(+.000000000000617809E0,+.000000000000156198E0), 00237 G4complex(-.000000000000043162E0,-.000000000000103385E0), 00238 G4complex(-.000000000000010133E0,+.000000000000018129E0), 00239 G4complex(+.000000000000003973E0,-.000000000000000709E0), 00240 G4complex(-.000000000000000632E0,-.000000000000000520E0), 00241 G4complex(+.000000000000000006E0,+.000000000000000172E0), 00242 G4complex(+.000000000000000027E0,-.000000000000000026E0), 00243 G4complex(-.000000000000000008E0,-.000000000000000000E0), 00244 G4complex(+.000000000000000001E0,+.000000000000000001E0)}; 00245 // ------------------------------------------------------------------------------------- 00246 G4double H=0.; // Prototype of the result 00247 if (ij) // I/K Bessel functions 00248 { 00249 if (ik) // I0/I1/EI0/EI1 Bessel functions (symmetric) 00250 { 00251 G4double V=std::abs(X); 00252 G4double CJ=0.; // Prototype of the element of the CI0/CI1 matrix 00253 if (V < 8.) 00254 { 00255 G4double HFV=HF*V; 00256 G4double Y=HFV*HFV; 00257 G4int V3=nu+1; 00258 G4int XL=V3+1; 00259 G4int XLI=XL+1; 00260 G4int XLD=XLI+1; 00261 G4int W1=XLD+1; 00262 G4double A0=1.; 00263 G4double DY=Y+Y; 00264 G4double A1=1.+DY/(XLI*V3); 00265 G4double A2=1.+Y*(4.+(DY+Y)/(XLD*XL))/(W1*V3); 00266 G4double B0=1.; 00267 G4double B1=1.-Y/XLI; 00268 G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1; 00269 G4int V1=3-XL; 00270 G4double V2=V3+V3; 00271 G4double C=0.; 00272 for (G4int N=3; N<=30; N++) 00273 { 00274 G4double C0=C; 00275 G4double FN=N; 00276 W1=W1+2; 00277 G4int W2=W1-1; 00278 G4int W3=W2-1; 00279 G4int W4=W3-1; 00280 G4int W5=W4-1; 00281 G4int W6=W5-1; 00282 V1=V1+1; 00283 V2=V2+1; 00284 V3=V3+1; 00285 G4double U1=FN*W4; 00286 G4double E=V3/(U1*W3); 00287 G4double U2=E*Y; 00288 G4double F1=1.+Y*V1/(U1*W1); 00289 G4double F2=(1.+Y*V2/(V3*W2*W5))*U2; 00290 G4double F3=-Y*Y*U2/(W4*W5*W5*W6); 00291 G4double A=F1*A2+F2*A1+F3*A0; 00292 G4double B=F1*B2+F2*B1+F3*B0; 00293 C=A/B; 00294 if (std::abs(C0-C) < EPS*std::abs(C)) break; 00295 A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B; 00296 } 00297 H=C; 00298 if (nu==1) H*=HF*X; 00299 if (ex) H*=std::exp(-V); 00300 } 00301 else 00302 { 00303 G4double P=16./V-1.; 00304 G4double ALFA=P+P; 00305 G4double B1=0.; 00306 G4double B2=0.; 00307 for (G4int I=npil; I>=0; I--) 00308 { 00309 if (!nu) CJ=CI0[I]; 00310 else CJ=CI1[I]; 00311 G4double B0=CJ+ALFA*B1-B2; 00312 B2=B1; 00313 B1=B0; 00314 } 00315 H=std::sqrt(RPI2/V)*(B1-P*B2); 00316 if (nu && X < 0.) H=-H; 00317 if (!ex) H*=std::exp(V); 00318 } 00319 } 00320 else // K0/K1/EK0/EK1 Bessel functions 00321 { 00322 #ifdef debug 00323 G4cout<<"G4BesIKJY: ------------>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; 00324 #endif 00325 G4double CK=0.; // Prototype of the element of the CI0/CI1 matrix 00326 if (X < 0.) 00327 { 00328 G4cout<<"G4BesIKJY::NegativeArg in K-BessFun X="<<X<<", n="<<nu<<",E="<<ex<<G4endl; 00329 return H; 00330 } 00331 else if (X < 1.) 00332 { 00333 #ifdef debug 00334 G4cout<<"G4BesIKJY: -->> [ X < 1 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; 00335 #endif 00336 G4double B=HF*X; 00337 G4double BK=-(std::log(B)+CE); 00338 G4double F=BK; 00339 G4double P=HF; 00340 G4double Q=HF; 00341 G4double C=1.; 00342 G4double D=B*B; 00343 G4double BK1=P; 00344 for (G4int N=1; N<=nat1; N++) // @@ "nat" can be increased 00345 { 00346 G4double FN=N; 00347 P/=FN; 00348 Q/=FN; 00349 F=(F+P+Q)/FN; 00350 C*=D/FN; 00351 G4double G=C*(P-FN*F); 00352 G4double R=C*F; 00353 BK=BK+R; 00354 BK1=BK1+G; 00355 if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break; 00356 } 00357 if (nu==1) H=BK1/B; 00358 else H=BK; 00359 if (ex) H*=std::exp(X); 00360 } 00361 else if (X < 5.) 00362 { 00363 #ifdef debug 00364 G4cout<<"G4BesIKJY: -->> [ X < 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; 00365 #endif 00366 G4int NUS=0; // @@ NU**2 for future NU>1 applications 00367 if (nu==1) NUS=1; 00368 G4double DNUS=NUS+NUS; 00369 G4double XN=DNUS+DNUS; 00370 G4double A=9.-XN; 00371 G4double B=25.-XN; 00372 G4double C=768*X*X; 00373 G4double HX=16*X; 00374 G4double C0=HX+HX+HX;; 00375 G4double A0=1.; 00376 G4double A1=(HX+7.+XN)/A; 00377 G4double A2=(C+C0*(XN+23.)+XN*(XN+62.)+129.)/(A*B); 00378 G4double B0=1.; 00379 G4double B1=(HX+9.-XN)/A; 00380 G4double B2=(C+C0*B)/(A*B)+1.; 00381 C=0.; 00382 for (G4int N=3; N<=nat2; N++) 00383 { 00384 C0=C; 00385 G4double FN=N; 00386 G4double FN2=FN+FN; 00387 G4double FNP=FN2+1.; 00388 G4double FN1=FN2-1.; 00389 G4double FNM=FN1-4.; 00390 G4double FN3=FN1/(FN2-3.); 00391 G4double FN4=12*FN*FN-(1.-XN); 00392 G4double FN5=16*FN1*X; 00393 G4double RAN=1./(FNP*FNP-XN); 00394 G4double F1=FN3*(FN4-20*FN)+FN5; 00395 G4double F2=28*FN-FN4-8.+FN5; 00396 G4double F3=FN3*(FNM*FNM-XN); 00397 A=(F1*A2+F2*A1+F3*A0)*RAN; 00398 B=(F1*B2+F2*B1+F3*B0)*RAN; 00399 C=A/B; 00400 if (std::abs(C0-C) < EPS*std::abs(C)) break; 00401 A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B; 00402 } 00403 H=C/std::sqrt(RPIH*X); 00404 if (!ex) H*=std::exp(-X); 00405 } 00406 else 00407 { 00408 #ifdef debug 00409 G4cout<<"G4BesIKJY: ->> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; 00410 #endif 00411 G4double P=10./X-1.; 00412 G4double ALFA=P+P; 00413 G4double B1=0.; 00414 G4double B2=0.; 00415 #ifdef debug 00416 G4cout<<"G4BesIKJY: ->> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; 00417 #endif 00418 for (G4int I=npkl; I>=0; I--) 00419 { 00420 if (!nu) CK=CK0[I]; 00421 else CK=CK1[I]; 00422 G4double B0=CK+ALFA*B1-B2; 00423 B2=B1; 00424 B1=B0; 00425 } 00426 H=std::sqrt(PIH/X)*(B1-P*B2); 00427 if (!ex) H*=std::exp(-X); 00428 } 00429 } 00430 } 00431 else 00432 { 00433 if (!ik && X < 0.) 00434 { 00435 G4cout<<"G4BesIKJY::NegativeArgument in Y BesselFunction X="<<X<<", nu="<<nu<<G4endl; 00436 return H; 00437 } 00438 G4double V=std::abs(X); 00439 if (!nu) // J0/Y0 Bessel functions 00440 { 00441 if (V < 8.) 00442 { 00443 G4double P=R32*V*V-1.; 00444 G4double ALFA=P+P; 00445 G4double B1=0.; 00446 G4double B2=0.; 00447 for (G4int IT=npkl; IT>=0; IT--) 00448 { 00449 G4double B0=CA[IT]+ALFA*B1-B2; 00450 B2=B1; 00451 B1=B0; 00452 } 00453 H=B1-P*B2; 00454 if (!ik) 00455 { 00456 B1=0.; 00457 B2=0.; 00458 for (G4int JT=npkl; JT>=0; JT--) 00459 { 00460 G4double B0=CB[JT]+ALFA*B1-B2; 00461 B2=B1; 00462 B1=B0; 00463 } 00464 H=RPIH*(CE+std::log(HF*X))*H+B1-P*B2; 00465 } 00466 } 00467 else 00468 { 00469 G4double P=10./V-1.; 00470 G4double ALFA=P+P; 00471 G4complex CB1(0.,0.); 00472 G4complex CB2(0.,0.); 00473 for (G4int IT=npjl; IT>=0; IT--) 00474 { 00475 G4complex CB0=CC[IT]+ALFA*CB1-CB2; 00476 CB2=CB1; 00477 CB1=CB0; 00478 } 00479 CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2); 00480 if (ik) H=real(CB1); 00481 else H=real(-CI*CB1); 00482 } 00483 } 00484 else // J1/Y1 Bessel functions 00485 { 00486 if (V < 8.) 00487 { 00488 G4double Y=R8*V; 00489 G4double Y2=Y*Y; 00490 G4double P=Y2+Y2-1.; 00491 G4double ALFA=P+P; 00492 G4double B1=0.; 00493 G4double B2=0.; 00494 for (G4int IT=npkl; IT>=0; IT--) 00495 { 00496 G4double B0=CD[IT]+ALFA*B1-B2; 00497 B2=B1; 00498 B1=B0; 00499 } 00500 H=Y*(B1-P*B2); 00501 if (!ik) 00502 { 00503 B1=0.; 00504 B2=0.; 00505 for (G4int JT=npkl; JT>=0; JT--) 00506 { 00507 G4double B0=EE[JT]+ALFA*B1-B2; 00508 B2=B1; 00509 B1=B0; 00510 } 00511 H=RPIH*((CE+std::log(HF*X))*H-1./X)+Y*(B1-B2); 00512 } 00513 } 00514 else 00515 { 00516 G4double P=10./V-1.; 00517 G4double ALFA=P+P; 00518 G4complex CB1(0.,0.); 00519 G4complex CB2(0.,0.); 00520 for (G4int IT=npjl; IT>=0; IT--) 00521 { 00522 G4complex CB0=CF[IT]+ALFA*CB1-CB2; 00523 CB2=CB1; 00524 CB1=CB0; 00525 } 00526 CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI3))*(CB1-P*CB2); 00527 if (ik) H=real(CB1); 00528 else H=real(-CI*CB1); 00529 } 00530 if (X < 0.) H=-H; 00531 } 00532 } 00533 return H; 00534 }
1.4.7