[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 | // |
---|
[1196] | 27 | // $Id: G4QBesIKJY.cc,v 1.4 2009/11/10 17:13:46 mkossov Exp $ |
---|
[1228] | 28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
[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 | // ------------------------------------------------------------------- |
---|
[1196] | 36 | |
---|
[819] | 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) |
---|
[1055] | 47 | { |
---|
| 48 | case BessI0: |
---|
[819] | 49 | nu=0; |
---|
| 50 | ik=true; |
---|
| 51 | ij=true; |
---|
[1055] | 52 | break; |
---|
| 53 | case BessI1: |
---|
[819] | 54 | nu=1; |
---|
| 55 | ik=true; |
---|
| 56 | ij=true; |
---|
[1055] | 57 | break; |
---|
| 58 | case EBessI0: |
---|
[819] | 59 | nu=0; |
---|
| 60 | ex=true; |
---|
| 61 | ik=true; |
---|
| 62 | ij=true; |
---|
[1055] | 63 | break; |
---|
| 64 | case EBessI1: |
---|
[819] | 65 | nu=1; |
---|
| 66 | ex=true; |
---|
| 67 | ik=true; |
---|
| 68 | ij=true; |
---|
[1055] | 69 | break; |
---|
| 70 | case BessJ0: |
---|
[819] | 71 | nu=0; |
---|
| 72 | ik=true; |
---|
| 73 | ij=false; |
---|
[1055] | 74 | break; |
---|
| 75 | case BessJ1: |
---|
[819] | 76 | nu=1; |
---|
| 77 | ik=true; |
---|
| 78 | ij=false; |
---|
[1055] | 79 | break; |
---|
| 80 | case BessK0: |
---|
[819] | 81 | nu=0; |
---|
| 82 | ik=false; |
---|
| 83 | ij=true; |
---|
[1055] | 84 | break; |
---|
| 85 | case BessK1: |
---|
[819] | 86 | nu=1; |
---|
| 87 | ik=false; |
---|
| 88 | ij=true; |
---|
[1055] | 89 | break; |
---|
| 90 | case EBessK0: |
---|
[819] | 91 | nu=0; |
---|
| 92 | ex=true; |
---|
| 93 | ik=false; |
---|
| 94 | ij=true; |
---|
[1055] | 95 | break; |
---|
| 96 | case EBessK1: |
---|
[819] | 97 | nu=1; |
---|
| 98 | ex=true; |
---|
| 99 | ik=false; |
---|
| 100 | ij=true; |
---|
[1055] | 101 | break; |
---|
| 102 | case BessY0: |
---|
[819] | 103 | nu=0; |
---|
| 104 | ik=false; |
---|
| 105 | ij=false; |
---|
[1055] | 106 | break; |
---|
| 107 | case BessY1: |
---|
[819] | 108 | nu=1; |
---|
| 109 | ik=false; |
---|
| 110 | ij=false; |
---|
[1055] | 111 | break; |
---|
[819] | 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.; |
---|
[1055] | 130 | static const G4double HF = Z1/2; |
---|
| 131 | static const G4double R8 = HF/4; |
---|
| 132 | static const G4double R32 = R8/4; |
---|
[819] | 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, |
---|
[1055] | 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}; |
---|
[819] | 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, |
---|
[1055] | 219 | -.022719244428417736E0,+.206644541017490520E0,-.086671697056948524E0, |
---|
| 220 | +.017636703003163134E0,-.002235619294485095E0,+.000197062302701541E0, |
---|
| 221 | -.000012885853299241E0,+.000000652847952359E0,-.000000026450737175E0, |
---|
| 222 | +.000000000878030117E0,-.000000000024343279E0,+.000000000000572612E0, |
---|
| 223 | -.000000000000011578E0,+.000000000000000203E0,-.000000000000000003E0}; |
---|
[819] | 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 |
---|
[1055] | 248 | if (ij) // I/K Bessel functions |
---|
[819] | 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); |
---|
[1055] | 267 | G4double B0=1.; |
---|
| 268 | G4double B1=1.-Y/XLI; |
---|
[819] | 269 | G4double B2=1.-Y*(1.-Y/(XLD+XLD))/W1; |
---|
| 270 | G4int V1=3-XL; |
---|
| 271 | G4double V2=V3+V3; |
---|
| 272 | G4double C=0.; |
---|
[1055] | 273 | for (G4int N=3; N<=30; N++) |
---|
[819] | 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; |
---|
[1055] | 283 | V1=V1+1; |
---|
[819] | 284 | V2=V2+1; |
---|
| 285 | V3=V3+1; |
---|
[1055] | 286 | G4double U1=FN*W4; |
---|
[819] | 287 | G4double E=V3/(U1*W3); |
---|
[1055] | 288 | G4double U2=E*Y; |
---|
[819] | 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; |
---|
[1055] | 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; |
---|
[819] | 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.; |
---|
[1055] | 305 | G4double ALFA=P+P; |
---|
| 306 | G4double B1=0.; |
---|
[819] | 307 | G4double B2=0.; |
---|
[1055] | 308 | for (G4int I=npil; I>=0; I--) |
---|
[819] | 309 | { |
---|
| 310 | if (!nu) CJ=CI0[I]; |
---|
| 311 | else CJ=CI1[I]; |
---|
[1055] | 312 | G4double B0=CJ+ALFA*B1-B2; |
---|
| 313 | B2=B1; |
---|
| 314 | B1=B0; |
---|
[819] | 315 | } |
---|
[1055] | 316 | H=std::sqrt(RPI2/V)*(B1-P*B2); |
---|
| 317 | if (nu && X < 0.) H=-H; |
---|
| 318 | if (!ex) H*=std::exp(V); |
---|
| 319 | } |
---|
[819] | 320 | } |
---|
| 321 | else // K0/K1/EK0/EK1 Bessel functions |
---|
[1055] | 322 | { |
---|
[819] | 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.; |
---|
[1055] | 343 | G4double D=B*B; |
---|
[819] | 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); |
---|
[1055] | 353 | G4double R=C*F; |
---|
[819] | 354 | BK=BK+R; |
---|
| 355 | BK1=BK1+G; |
---|
[1055] | 356 | if (BK1*R+std::abs(G)*BK < EPS*BK*BK1) break; |
---|
[819] | 357 | } |
---|
| 358 | if (nu==1) H=BK1/B; |
---|
| 359 | else H=BK; |
---|
| 360 | if (ex) H*=std::exp(X); |
---|
| 361 | } |
---|
[1055] | 362 | else if (X < 5.) |
---|
[819] | 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; |
---|
[1055] | 372 | G4double B=25.-XN; |
---|
[819] | 373 | G4double C=768*X*X; |
---|
[1055] | 374 | G4double HX=16*X; |
---|
| 375 | G4double C0=HX+HX+HX;; |
---|
[819] | 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.; |
---|
[1055] | 382 | C=0.; |
---|
| 383 | for (G4int N=3; N<=nat2; N++) |
---|
[819] | 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); |
---|
[1055] | 393 | G4double FN5=16*FN1*X; |
---|
[819] | 394 | G4double RAN=1./(FNP*FNP-XN); |
---|
[1055] | 395 | G4double F1=FN3*(FN4-20*FN)+FN5; |
---|
[819] | 396 | G4double F2=28*FN-FN4-8.+FN5; |
---|
| 397 | G4double F3=FN3*(FNM*FNM-XN); |
---|
[1055] | 398 | A=(F1*A2+F2*A1+F3*A0)*RAN; |
---|
[819] | 399 | B=(F1*B2+F2*B1+F3*B0)*RAN; |
---|
| 400 | C=A/B; |
---|
[1055] | 401 | if (std::abs(C0-C) < EPS*std::abs(C)) break; |
---|
| 402 | A0=A1; A1=A2; A2=A; B0=B1; B1=B2; B2=B; |
---|
[819] | 403 | } |
---|
[1055] | 404 | H=C/std::sqrt(RPIH*X); |
---|
[819] | 405 | if (!ex) H*=std::exp(-X); |
---|
[1055] | 406 | } |
---|
[819] | 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; |
---|
[1055] | 414 | G4double B1=0.; |
---|
| 415 | G4double B2=0.; |
---|
[819] | 416 | #ifdef debug |
---|
| 417 | G4cout<<"G4BesIKJY: >>> [ X >= 5 ] is called, X="<<X<<",n="<<nu<<",E="<<ex<<G4endl; |
---|
| 418 | #endif |
---|
[1055] | 419 | for (G4int I=npkl; I>=0; I--) |
---|
[819] | 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; |
---|
[1055] | 426 | } |
---|
[819] | 427 | H=std::sqrt(PIH/X)*(B1-P*B2); |
---|
| 428 | if (!ex) H*=std::exp(-X); |
---|
| 429 | } |
---|
[1055] | 430 | } |
---|
[819] | 431 | } |
---|
| 432 | else |
---|
[1055] | 433 | { |
---|
[819] | 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.; |
---|
[1055] | 459 | for (G4int JT=npkl; JT>=0; JT--) |
---|
[819] | 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 | } |
---|
[1055] | 467 | } |
---|
[819] | 468 | else |
---|
[1055] | 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--) |
---|
[819] | 475 | { |
---|
| 476 | G4complex CB0=CC[IT]+ALFA*CB1-CB2; |
---|
| 477 | CB2=CB1; |
---|
| 478 | CB1=CB0; |
---|
[1055] | 479 | } |
---|
[819] | 480 | CB1=std::sqrt(RPIH/V)*std::exp(CI*(V-PI4))*(CB1-P*CB2); |
---|
[1055] | 481 | if (ik) H=real(CB1); |
---|
[819] | 482 | else H=real(-CI*CB1); |
---|
[1055] | 483 | } |
---|
[819] | 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; |
---|
[1055] | 498 | B2=B1; |
---|
[819] | 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 | } |
---|
[1055] | 514 | } |
---|
[819] | 515 | else |
---|
[1055] | 516 | { |
---|
[819] | 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; |
---|
[1055] | 524 | CB2=CB1; |
---|
[819] | 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); |
---|
[1055] | 530 | } |
---|
[819] | 531 | if (X < 0.) H=-H; |
---|
| 532 | } |
---|
| 533 | } |
---|
| 534 | return H; |
---|
| 535 | } |
---|