G4QBesIKJY.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //      ---------------- G4QBesIKJY ----------------
00030 //             by Mikhail Kossov, Sept 1999.
00031 //      class  for Bessel I0/I1 and K0/K1 functions in CHIPS Model
00032 // -------------------------------------------------------------------
00033 // Short description: Bessel functions class (can be substituted)
00034 // -------------------------------------------------------------------
00035 
00036 //#define debug
00037 //#define pdebug
00038 
00039 #include "G4QBesIKJY.hh"
00040 
00041 // Constructor
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;            // 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 }

Generated on Mon May 27 17:49:30 2013 for Geant4 by  doxygen 1.4.7