[819] | 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 | // |
---|
[1055] | 27 | // $Id: G4QBesIKJY.cc,v 1.3 2009/02/23 09:49:24 mkossov Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
[819] | 29 | // |
---|
| 30 | // ---------------- G4QBesIKJY ---------------- |
---|
| 31 | // by Mikhail Kossov, Sept 1999. |
---|
| 32 | // class for Bessel I0/I1 and K0/K1 functions in CHIPS Model |
---|
| 33 | // ------------------------------------------------------------------- |
---|
[1055] | 34 | // Short description: Bessel functions class (can be substituted) |
---|
| 35 | // ------------------------------------------------------------------- |
---|
[819] | 36 | //#define debug |
---|
| 37 | //#define pdebug |
---|
| 38 | |
---|
| 39 | #include "G4QBesIKJY.hh" |
---|
| 40 | |
---|
| 41 | // Constructor |
---|
| 42 | G4QBesIKJY::G4QBesIKJY(G4QBIType type) |
---|
| 43 | { |
---|
| 44 | ex=false; |
---|
| 45 | switch (type) |
---|
[1055] | 46 | { |
---|
| 47 | case BessI0: |
---|
[819] | 48 | nu=0; |
---|
| 49 | ik=true; |
---|
| 50 | ij=true; |
---|
[1055] | 51 | break; |
---|
| 52 | case BessI1: |
---|
[819] | 53 | nu=1; |
---|
| 54 | ik=true; |
---|
| 55 | ij=true; |
---|
[1055] | 56 | break; |
---|
| 57 | case EBessI0: |
---|
[819] | 58 | nu=0; |
---|
| 59 | ex=true; |
---|
| 60 | ik=true; |
---|
| 61 | ij=true; |
---|
[1055] | 62 | break; |
---|
| 63 | case EBessI1: |
---|
[819] | 64 | nu=1; |
---|
| 65 | ex=true; |
---|
| 66 | ik=true; |
---|
| 67 | ij=true; |
---|
[1055] | 68 | break; |
---|
| 69 | case BessJ0: |
---|
[819] | 70 | nu=0; |
---|
| 71 | ik=true; |
---|
| 72 | ij=false; |
---|
[1055] | 73 | break; |
---|
| 74 | case BessJ1: |
---|
[819] | 75 | nu=1; |
---|
| 76 | ik=true; |
---|
| 77 | ij=false; |
---|
[1055] | 78 | break; |
---|
| 79 | case BessK0: |
---|
[819] | 80 | nu=0; |
---|
| 81 | ik=false; |
---|
| 82 | ij=true; |
---|
[1055] | 83 | break; |
---|
| 84 | case BessK1: |
---|
[819] | 85 | nu=1; |
---|
| 86 | ik=false; |
---|
| 87 | ij=true; |
---|
[1055] | 88 | break; |
---|
| 89 | case EBessK0: |
---|
[819] | 90 | nu=0; |
---|
| 91 | ex=true; |
---|
| 92 | ik=false; |
---|
| 93 | ij=true; |
---|
[1055] | 94 | break; |
---|
| 95 | case EBessK1: |
---|
[819] | 96 | nu=1; |
---|
| 97 | ex=true; |
---|
| 98 | ik=false; |
---|
| 99 | ij=true; |
---|
[1055] | 100 | break; |
---|
| 101 | case BessY0: |
---|
[819] | 102 | nu=0; |
---|
| 103 | ik=false; |
---|
| 104 | ij=false; |
---|
[1055] | 105 | break; |
---|
| 106 | case BessY1: |
---|
[819] | 107 | nu=1; |
---|
| 108 | ik=false; |
---|
| 109 | ij=false; |
---|
[1055] | 110 | break; |
---|
[819] | 111 | } |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | G4QBesIKJY::~G4QBesIKJY(){;} |
---|
| 115 | |
---|
| 116 | G4double 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.; |
---|
[1055] | 129 | static const G4double HF = Z1/2; |
---|
| 130 | static const G4double R8 = HF/4; |
---|
| 131 | static const G4double R32 = R8/4; |
---|
[819] | 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, |
---|
[1055] | 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}; |
---|
[819] | 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, |
---|
[1055] | 218 | -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0, |
---|
| 219 | +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0, |
---|
| 220 | -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0, |
---|
| 221 | +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0, |
---|
| 222 | -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0}; |
---|
[819] | 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 |
---|
[1055] | 247 | if (ij) // I/K Bessel functions |
---|
[819] | 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); |
---|
[1055] | 266 | G4double B0=1.; |
---|
| 267 | G4double B1=1.-Y/XLI; |
---|
[819] | 268 | G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1; |
---|
| 269 | G4int V1=3-XL; |
---|
| 270 | G4double V2=V3+V3; |
---|
| 271 | G4double C=0.; |
---|
[1055] | 272 | for (G4int N=3; N<=30; N++) |
---|
[819] | 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; |
---|
[1055] | 282 | V1=V1+1; |
---|
[819] | 283 | V2=V2+1; |
---|
| 284 | V3=V3+1; |
---|
[1055] | 285 | G4double U1=FN*W4; |
---|
[819] | 286 | G4double E=V3/(U1*W3); |
---|
[1055] | 287 | G4double U2=E*Y; |
---|
[819] | 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; |
---|
[1055] | 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; |
---|
[819] | 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.; |
---|
[1055] | 304 | G4double ALFA=P+P; |
---|
| 305 | G4double B1=0.; |
---|
[819] | 306 | G4double B2=0.; |
---|
[1055] | 307 | for (G4int I=npil; I>=0; I--) |
---|
[819] | 308 | { |
---|
| 309 | if (!nu) CJ=CI0[I]; |
---|
| 310 | else CJ=CI1[I]; |
---|
[1055] | 311 | G4double B0=CJ+ALFA*B1-B2; |
---|
| 312 | B2=B1; |
---|
| 313 | B1=B0; |
---|
[819] | 314 | } |
---|
[1055] | 315 | H=std::sqrt(RPI2/V)*(B1-P*B2); |
---|
| 316 | if (nu && X < 0.) H=-H; |
---|
| 317 | if (!ex) H*=std::exp(V); |
---|
| 318 | } |
---|
[819] | 319 | } |
---|
| 320 | else // K0/K1/EK0/EK1 Bessel functions |
---|
[1055] | 321 | { |
---|
[819] | 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.; |
---|
[1055] | 342 | G4double D=B*B; |
---|
[819] | 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); |
---|
[1055] | 352 | G4double R=C*F; |
---|
[819] | 353 | BK=BK+R; |
---|
| 354 | BK1=BK1+G; |
---|
[1055] | 355 | if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break; |
---|
[819] | 356 | } |
---|
| 357 | if (nu==1) H=BK1/B; |
---|
| 358 | else H=BK; |
---|
| 359 | if (ex) H*=std::exp(X); |
---|
| 360 | } |
---|
[1055] | 361 | else if (X < 5.) |
---|
[819] | 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; |
---|
[1055] | 371 | G4double B=25.-XN; |
---|
[819] | 372 | G4double C=768*X*X; |
---|
[1055] | 373 | G4double HX=16*X; |
---|
| 374 | G4double C0=HX+HX+HX;; |
---|
[819] | 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.; |
---|
[1055] | 381 | C=0.; |
---|
| 382 | for (G4int N=3; N<=nat2; N++) |
---|
[819] | 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); |
---|
[1055] | 392 | G4double FN5=16*FN1*X; |
---|
[819] | 393 | G4double RAN=1./(FNP*FNP-XN); |
---|
[1055] | 394 | G4double F1=FN3*(FN4-20*FN)+FN5; |
---|
[819] | 395 | G4double F2=28*FN-FN4-8.+FN5; |
---|
| 396 | G4double F3=FN3*(FNM*FNM-XN); |
---|
[1055] | 397 | A=(F1*A2+F2*A1+F3*A0)*RAN; |
---|
[819] | 398 | B=(F1*B2+F2*B1+F3*B0)*RAN; |
---|
| 399 | C=A/B; |
---|
[1055] | 400 | if (std::abs(C0-C) < EPS*std::abs(C)) break; |
---|
| 401 | A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B; |
---|
[819] | 402 | } |
---|
[1055] | 403 | H=C/std::sqrt(RPIH*X); |
---|
[819] | 404 | if (!ex) H*=std::exp(-X); |
---|
[1055] | 405 | } |
---|
[819] | 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; |
---|
[1055] | 413 | G4double B1=0.; |
---|
| 414 | G4double B2=0.; |
---|
[819] | 415 | #ifdef debug |
---|
| 416 | G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; |
---|
| 417 | #endif |
---|
[1055] | 418 | for (G4int I=npkl; I>=0; I--) |
---|
[819] | 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; |
---|
[1055] | 425 | } |
---|
[819] | 426 | H=std::sqrt(PIH/X)*(B1-P*B2); |
---|
| 427 | if (!ex) H*=std::exp(-X); |
---|
| 428 | } |
---|
[1055] | 429 | } |
---|
[819] | 430 | } |
---|
| 431 | else |
---|
[1055] | 432 | { |
---|
[819] | 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.; |
---|
[1055] | 458 | for (G4int JT=npkl; JT>=0; JT--) |
---|
[819] | 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 | } |
---|
[1055] | 466 | } |
---|
[819] | 467 | else |
---|
[1055] | 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--) |
---|
[819] | 474 | { |
---|
| 475 | G4complex CB0=CC[IT]+ALFA*CB1-CB2; |
---|
| 476 | CB2=CB1; |
---|
| 477 | CB1=CB0; |
---|
[1055] | 478 | } |
---|
[819] | 479 | CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2); |
---|
[1055] | 480 | if (ik) H=real(CB1); |
---|
[819] | 481 | else H=real(-CI*CB1); |
---|
[1055] | 482 | } |
---|
[819] | 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; |
---|
[1055] | 497 | B2=B1; |
---|
[819] | 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 | } |
---|
[1055] | 513 | } |
---|
[819] | 514 | else |
---|
[1055] | 515 | { |
---|
[819] | 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; |
---|
[1055] | 523 | CB2=CB1; |
---|
[819] | 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); |
---|
[1055] | 529 | } |
---|
[819] | 530 | if (X < 0.) H=-H; |
---|
| 531 | } |
---|
| 532 | } |
---|
| 533 | return H; |
---|
| 534 | } |
---|