00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4QBesIKJY.hh"
00040
00041
00042 G4QBesIKJY::G4QBesIKJY(G4QBIType type)
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 }
00113
00114 G4QBesIKJY::~G4QBesIKJY(){;}
00115
00116 G4double G4QBesIKJY::operator() (G4double X) const
00117 {
00118 static const G4int nat1 = 15;
00119 static const G4int nat2 = nat1+nat1;
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;
00136 static const G4double PI3 = PIH+PI4;
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.;
00247 if (ij)
00248 {
00249 if (ik)
00250 {
00251 G4double V=std::abs(X);
00252 G4double CJ=0.;
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
00321 {
00322 #ifdef debug
00323 G4cout<<"G4BesIKJY: ------------>> K is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl;
00324 #endif
00325 G4double CK=0.;
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++)
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;
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)
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
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 }