Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * 20 // * * 21 // * Parts of this code which have been devel 21 // * Parts of this code which have been developed by QinetiQ Ltd * 22 // * under contract to the European Space Agen 22 // * under contract to the European Space Agency (ESA) are the * 23 // * intellectual property of ESA. Rights to u 23 // * intellectual property of ESA. Rights to use, copy, modify and * 24 // * redistribute this software for general pu 24 // * redistribute this software for general public use are granted * 25 // * in compliance with any licensing, distrib 25 // * in compliance with any licensing, distribution and development * 26 // * policy adopted by the Geant4 Collaboratio 26 // * policy adopted by the Geant4 Collaboration. This code has been * 27 // * written by QinetiQ Ltd for the European S 27 // * written by QinetiQ Ltd for the European Space Agency, under ESA * 28 // * contract 17191/03/NL/LvH (Aurora Programm 28 // * contract 17191/03/NL/LvH (Aurora Programme). * 29 // * 29 // * * 30 // * By using, copying, modifying or distri 30 // * By using, copying, modifying or distributing the software (or * 31 // * any work based on the software) you ag 31 // * any work based on the software) you agree to acknowledge its * 32 // * use in resulting scientific publicati 32 // * use in resulting scientific publications, and indicate your * 33 // * acceptance of all terms of the Geant4 Sof 33 // * acceptance of all terms of the Geant4 Software license. * 34 // ******************************************* 34 // ******************************************************************** 35 // 35 // 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37 // 37 // 38 // MODULE: G4Bessel.cc 38 // MODULE: G4Bessel.cc 39 // 39 // 40 // Version: B.1 40 // Version: B.1 41 // Date: 15/04/04 41 // Date: 15/04/04 42 // Author: P R Truscott 42 // Author: P R Truscott 43 // Organisation: QinetiQ Ltd, UK 43 // Organisation: QinetiQ Ltd, UK 44 // Customer: ESA/ESTEC, NOORDWIJK 44 // Customer: ESA/ESTEC, NOORDWIJK 45 // Contract: 17191/03/NL/LvH 45 // Contract: 17191/03/NL/LvH 46 // 46 // 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 48 // 48 // 49 // CHANGE HISTORY 49 // CHANGE HISTORY 50 // -------------- 50 // -------------- 51 // 51 // 52 // 18 Noevmber 2003, P R Truscott, QinetiQ Ltd 52 // 18 Noevmber 2003, P R Truscott, QinetiQ Ltd, UK 53 // Created. 53 // Created. 54 // 54 // 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, U 55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK 56 // Beta release 56 // Beta release 57 // 57 // 58 // 06 August 2015, A. Ribon, CERN << 59 // Migrated to G4Exp, G4Log and G4Pow. << 60 // << 61 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 58 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 62 ////////////////////////////////////////////// 59 //////////////////////////////////////////////////////////////////////////////// 63 // 60 // 64 #include "G4Bessel.hh" 61 #include "G4Bessel.hh" 65 #include "G4PhysicalConstants.hh" << 66 << 67 #include "G4Exp.hh" << 68 #include "G4Log.hh" << 69 #include "G4Pow.hh" << 70 ////////////////////////////////////////////// 62 //////////////////////////////////////////////////////////////////////////////// 71 // 63 // 72 G4Bessel::G4Bessel () 64 G4Bessel::G4Bessel () 73 {;} 65 {;} 74 ////////////////////////////////////////////// 66 //////////////////////////////////////////////////////////////////////////////// 75 // 67 // 76 G4Bessel::~G4Bessel () 68 G4Bessel::~G4Bessel () 77 {;} 69 {;} 78 ////////////////////////////////////////////// 70 //////////////////////////////////////////////////////////////////////////////// 79 // 71 // 80 G4double G4Bessel::I0 (G4double x) 72 G4double G4Bessel::I0 (G4double x) 81 { 73 { 82 const G4double P1 = 1.0, 74 const G4double P1 = 1.0, 83 P2 = 3.5156229, 75 P2 = 3.5156229, 84 P3 = 3.0899424, 76 P3 = 3.0899424, 85 P4 = 1.2067492, 77 P4 = 1.2067492, 86 P5 = 0.2659732, 78 P5 = 0.2659732, 87 P6 = 0.0360768, 79 P6 = 0.0360768, 88 P7 = 0.0045813; 80 P7 = 0.0045813; 89 const G4double Q1 = 0.39894228, 81 const G4double Q1 = 0.39894228, 90 Q2 = 0.01328592, 82 Q2 = 0.01328592, 91 Q3 = 0.00225319, 83 Q3 = 0.00225319, 92 Q4 =-0.00157565, 84 Q4 =-0.00157565, 93 Q5 = 0.00916281, 85 Q5 = 0.00916281, 94 Q6 =-0.02057706, 86 Q6 =-0.02057706, 95 Q7 = 0.02635537, 87 Q7 = 0.02635537, 96 Q8 =-0.01647633, 88 Q8 =-0.01647633, 97 Q9 = 0.00392377; 89 Q9 = 0.00392377; 98 90 99 G4double I = 0.0; 91 G4double I = 0.0; 100 if (std::abs(x) < 3.75) << 92 if (std::fabs(x) < 3.75) 101 { 93 { 102 G4double y = G4Pow::GetInstance()->powN(x/ << 94 G4double y = std::pow(x/3.75, 2.0); 103 I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7) 95 I = P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))); 104 } 96 } 105 else 97 else 106 { 98 { 107 G4double ax = std::abs(x); << 99 G4double ax = std::fabs(x); 108 G4double y = 3.75/ax; 100 G4double y = 3.75/ax; 109 I = G4Exp(ax) / std::sqrt(ax) * << 101 I = std::exp(ax) / std::sqrt(ax) * 110 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+ 102 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9)))))))); 111 } 103 } 112 return I; 104 return I; 113 } 105 } 114 ////////////////////////////////////////////// 106 //////////////////////////////////////////////////////////////////////////////// 115 // 107 // 116 G4double G4Bessel::K0 (G4double x) 108 G4double G4Bessel::K0 (G4double x) 117 { 109 { 118 const G4double P1 =-0.57721566, 110 const G4double P1 =-0.57721566, 119 P2 = 0.42278420, 111 P2 = 0.42278420, 120 P3 = 0.23069756, 112 P3 = 0.23069756, 121 P4 = 0.03488590, 113 P4 = 0.03488590, 122 P5 = 0.00262698, 114 P5 = 0.00262698, 123 P6 = 0.00010750, 115 P6 = 0.00010750, 124 P7 = 0.00000740; 116 P7 = 0.00000740; 125 const G4double Q1 = 1.25331414, 117 const G4double Q1 = 1.25331414, 126 Q2 =-0.07832358, 118 Q2 =-0.07832358, 127 Q3 = 0.02189568, 119 Q3 = 0.02189568, 128 Q4 =-0.01062446, 120 Q4 =-0.01062446, 129 Q5 = 0.00587872, 121 Q5 = 0.00587872, 130 Q6 =-0.00251540, 122 Q6 =-0.00251540, 131 Q7 = 0.00053208; 123 Q7 = 0.00053208; 132 124 133 G4double K = 0.0; 125 G4double K = 0.0; 134 if (x <= 2.0) 126 if (x <= 2.0) 135 { 127 { 136 G4double y = x * x / 4.0; 128 G4double y = x * x / 4.0; 137 K = (-G4Log(x/2.0)) * I0(x) + << 129 K = (-std::log(x/2.0)) * I0(x) + 138 P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))) 130 P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7))))); 139 } 131 } 140 else 132 else 141 { 133 { 142 G4double y = 2.0 / x; 134 G4double y = 2.0 / x; 143 K = G4Exp(-x) / std::sqrt(x) * << 135 K = std::exp(-x) / std::sqrt(x) * 144 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7)) 136 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7)))))); 145 } 137 } 146 return K; 138 return K; 147 } 139 } 148 ////////////////////////////////////////////// 140 //////////////////////////////////////////////////////////////////////////////// 149 // 141 // 150 G4double G4Bessel::I1 (G4double x) 142 G4double G4Bessel::I1 (G4double x) 151 { 143 { 152 const G4double P1 = 0.5, 144 const G4double P1 = 0.5, 153 P2 = 0.87890594, 145 P2 = 0.87890594, 154 P3 = 0.51498869, 146 P3 = 0.51498869, 155 P4 = 0.15084934, 147 P4 = 0.15084934, 156 P5 = 0.02658733, 148 P5 = 0.02658733, 157 P6 = 0.00301532, 149 P6 = 0.00301532, 158 P7 = 0.00032411; 150 P7 = 0.00032411; 159 const G4double Q1 = 0.39894228, 151 const G4double Q1 = 0.39894228, 160 Q2 =-0.03988024, 152 Q2 =-0.03988024, 161 Q3 =-0.00362018, 153 Q3 =-0.00362018, 162 Q4 = 0.00163801, 154 Q4 = 0.00163801, 163 Q5 =-0.01031555, 155 Q5 =-0.01031555, 164 Q6 = 0.02282967, 156 Q6 = 0.02282967, 165 Q7 =-0.02895312, 157 Q7 =-0.02895312, 166 Q8 = 0.01787654, 158 Q8 = 0.01787654, 167 Q9 =-0.00420059; 159 Q9 =-0.00420059; 168 160 169 G4double I = 0.0; 161 G4double I = 0.0; 170 if (std::abs(x) < 3.75) << 162 if (std::fabs(x) < 3.75) 171 { 163 { 172 G4double ax = std::abs(x); << 164 G4double ax = std::fabs(x); 173 G4double y = G4Pow::GetInstance()->powN(x/ << 165 G4double y = std::pow(x/3.75, 2.0); 174 I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y 166 I = ax*(P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))))); 175 } 167 } 176 else 168 else 177 { 169 { 178 G4double ax = std::abs(x); << 170 G4double ax = std::fabs(x); 179 G4double y = 3.75/ax; 171 G4double y = 3.75/ax; 180 I = G4Exp(ax) / std::sqrt(ax) * << 172 I = std::exp(ax) / std::sqrt(ax) * 181 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+ 173 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*(Q7+y*(Q8+y*Q9)))))))); 182 } 174 } 183 if (x < 0.0) I = -I; 175 if (x < 0.0) I = -I; 184 return I; 176 return I; 185 } 177 } 186 ////////////////////////////////////////////// 178 //////////////////////////////////////////////////////////////////////////////// 187 // 179 // 188 G4double G4Bessel::K1 (G4double x) 180 G4double G4Bessel::K1 (G4double x) 189 { 181 { 190 const G4double P1 = 1.0, 182 const G4double P1 = 1.0, 191 P2 = 0.15443144, 183 P2 = 0.15443144, 192 P3 =-0.67278579, 184 P3 =-0.67278579, 193 P4 =-0.18156897, 185 P4 =-0.18156897, 194 P5 =-0.01919402, 186 P5 =-0.01919402, 195 P6 =-0.00110404, 187 P6 =-0.00110404, 196 P7 =-0.00004686; 188 P7 =-0.00004686; 197 const G4double Q1 = 1.25331414, 189 const G4double Q1 = 1.25331414, 198 Q2 = 0.23498619, 190 Q2 = 0.23498619, 199 Q3 =-0.03655620, 191 Q3 =-0.03655620, 200 Q4 = 0.01504268, 192 Q4 = 0.01504268, 201 Q5 =-0.00780353, 193 Q5 =-0.00780353, 202 Q6 = 0.00325614, 194 Q6 = 0.00325614, 203 Q7 =-0.00068245; 195 Q7 =-0.00068245; 204 196 205 G4double K = 0.0; 197 G4double K = 0.0; 206 if (x <= 2.0) 198 if (x <= 2.0) 207 { 199 { 208 G4double y = x * x / 4.0; 200 G4double y = x * x / 4.0; 209 K = G4Log(x/2.0)*I1(x) + 1.0/x * << 201 K = std::log(x/2.0)*I1(x) + 1.0/x * 210 (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)) 202 (P1+y*(P2+y*(P3+y*(P4+y*(P5+y*(P6+y*P7)))))); 211 } 203 } 212 else 204 else 213 { 205 { 214 G4double y = 2.0 / x; 206 G4double y = 2.0 / x; 215 K = G4Exp(-x) / std::sqrt(x) * << 207 K = std::exp(-x) / std::sqrt(x) * 216 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7)) 208 (Q1+y*(Q2+y*(Q3+y*(Q4+y*(Q5+y*(Q6+y*Q7)))))); 217 } 209 } 218 return K; 210 return K; 219 } 211 } 220 ////////////////////////////////////////////// 212 //////////////////////////////////////////////////////////////////////////////// 221 // 213 // 222 G4double G4Bessel::pI0 (G4double x) 214 G4double G4Bessel::pI0 (G4double x) 223 { 215 { 224 const G4double A0 = 0.1250000000000E+00, 216 const G4double A0 = 0.1250000000000E+00, 225 A1 = 7.0312500000000E-02, 217 A1 = 7.0312500000000E-02, 226 A2 = 7.3242187500000E-02, 218 A2 = 7.3242187500000E-02, 227 A3 = 1.1215209960938E-01, 219 A3 = 1.1215209960938E-01, 228 A4 = 2.2710800170898E-01, 220 A4 = 2.2710800170898E-01, 229 A5 = 5.7250142097473E-01, 221 A5 = 5.7250142097473E-01, 230 A6 = 1.7277275025845E+00, 222 A6 = 1.7277275025845E+00, 231 A7 = 6.0740420012735E+00, 223 A7 = 6.0740420012735E+00, 232 A8 = 2.4380529699556E+01, 224 A8 = 2.4380529699556E+01, 233 A9 = 1.1001714026925E+02, 225 A9 = 1.1001714026925E+02, 234 A10 = 5.5133589612202E+02, 226 A10 = 5.5133589612202E+02, 235 A11 = 3.0380905109224E+03; 227 A11 = 3.0380905109224E+03; 236 228 237 G4double I = 0.0; 229 G4double I = 0.0; 238 if (x == 0.0) 230 if (x == 0.0) 239 { 231 { 240 I = 1.0; 232 I = 1.0; 241 } 233 } 242 else if (x < 18.0) 234 else if (x < 18.0) 243 { 235 { 244 I = 1.0; 236 I = 1.0; 245 G4double y = x * x; 237 G4double y = x * x; 246 G4double q = 1.0; << 238 G4double s = 1.0; 247 for (G4int i=1; i<101; i++) 239 for (G4int i=1; i<101; i++) 248 { 240 { 249 q *= 0.25 * y / i / i; << 241 s *= 0.25 * y / i / i; 250 I += q; << 242 I += s; 251 if (std::abs(q/I) < 1.0E-15) break; << 243 if (std::abs(s/I) < 1.0E-15) break; 252 } 244 } 253 } 245 } 254 else 246 else 255 { 247 { 256 G4double y = 1.0 / x; 248 G4double y = 1.0 / x; 257 I = G4Exp(x) / std::sqrt(twopi*x) * << 249 I = std::exp(x) / std::sqrt(twopi*x) * 258 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*( 250 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11)))))))))))); 259 } 251 } 260 252 261 return I; 253 return I; 262 } 254 } 263 ////////////////////////////////////////////// 255 //////////////////////////////////////////////////////////////////////////////// 264 // 256 // 265 G4double G4Bessel::pI1 (G4double x) 257 G4double G4Bessel::pI1 (G4double x) 266 { 258 { 267 const G4double A0 = -0.3750000000000E+00, 259 const G4double A0 = -0.3750000000000E+00, 268 A1 = -1.1718750000000E-01, 260 A1 = -1.1718750000000E-01, 269 A2 = -1.0253906250000E-01, 261 A2 = -1.0253906250000E-01, 270 A3 = -1.4419555664063E-01, 262 A3 = -1.4419555664063E-01, 271 A4 = -2.775764465332E-01, 263 A4 = -2.775764465332E-01, 272 A5 = -6.7659258842468E-01, 264 A5 = -6.7659258842468E-01, 273 A6 = -1.9935317337513E+00, 265 A6 = -1.9935317337513E+00, 274 A7 = -6.8839142681099E+00, 266 A7 = -6.8839142681099E+00, 275 A8 = -2.7248827311269E+01, 267 A8 = -2.7248827311269E+01, 276 A9 = -1.2159789187654E+02, 268 A9 = -1.2159789187654E+02, 277 A10 = -6.0384407670507E+02, 269 A10 = -6.0384407670507E+02, 278 A11 = -3.3022722944809E+03; 270 A11 = -3.3022722944809E+03; 279 271 280 G4double I = 0.0; 272 G4double I = 0.0; 281 if (x == 0.0) 273 if (x == 0.0) 282 { 274 { 283 I = 0.0; 275 I = 0.0; 284 } 276 } 285 else if (x < 18.0) 277 else if (x < 18.0) 286 { 278 { 287 I = 1.0; 279 I = 1.0; 288 G4double y = x * x; 280 G4double y = x * x; 289 G4double q = 1.0; << 281 G4double s = 1.0; 290 for (G4int i=1; i<101; i++) 282 for (G4int i=1; i<101; i++) 291 { 283 { 292 q *= 0.25 * y / i / (i+1.0); << 284 s *= 0.25 * y / i / (i+1.0); 293 I += q; << 285 I += s; 294 if (std::abs(q/I) < 1.0E-15) break; << 286 if (std::abs(s/I) < 1.0E-15) break; 295 } 287 } 296 I *= 0.5 * x; 288 I *= 0.5 * x; 297 289 298 } 290 } 299 else 291 else 300 { 292 { 301 G4double y = 1.0 / x; 293 G4double y = 1.0 / x; 302 I = G4Exp(x) / std::sqrt(twopi*x) * << 294 I = std::exp(x) / std::sqrt(twopi*x) * 303 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*( 295 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*(A7+y*(A8+y*(A9+y*(A10+y*A11)))))))))))); 304 } 296 } 305 297 306 return I; 298 return I; 307 } 299 } 308 ////////////////////////////////////////////// 300 //////////////////////////////////////////////////////////////////////////////// 309 // 301 // 310 G4double G4Bessel::pK0 (G4double x) 302 G4double G4Bessel::pK0 (G4double x) 311 { 303 { 312 const G4double A0 = 0.1250000000000E+00, 304 const G4double A0 = 0.1250000000000E+00, 313 A1 = 0.2109375000000E+00, 305 A1 = 0.2109375000000E+00, 314 A2 = 1.0986328125000E+00, 306 A2 = 1.0986328125000E+00, 315 A3 = 1.1775970458984E+01, 307 A3 = 1.1775970458984E+01, 316 A4 = 2.1461706161499E+02, 308 A4 = 2.1461706161499E+02, 317 A5 = 5.9511522710323E+03, 309 A5 = 5.9511522710323E+03, 318 A6 = 2.3347645606175E+05, 310 A6 = 2.3347645606175E+05, 319 A7 = 1.2312234987631E+07; 311 A7 = 1.2312234987631E+07; 320 312 321 G4double K = 0.0; 313 G4double K = 0.0; 322 if (x == 0.0) 314 if (x == 0.0) 323 { 315 { 324 K = 1.0E+307; 316 K = 1.0E+307; 325 } 317 } 326 else if (x < 9.0) 318 else if (x < 9.0) 327 { 319 { 328 G4double y = x * x; 320 G4double y = x * x; 329 G4double C = -G4Log(x/2.0) - 0.57721566490 << 321 G4double C = -std::log(x/2.0) - 0.5772156649015329; 330 G4double q = 1.0; << 322 G4double s = 1.0; 331 G4double t = 0.0; 323 G4double t = 0.0; 332 for (G4int i=1; i<51; i++) 324 for (G4int i=1; i<51; i++) 333 { 325 { 334 q *= 0.25 * y / i / i; << 326 s *= 0.25 * y / i / i; 335 t += 1.0 / i ; 327 t += 1.0 / i ; 336 K += q * (t+C); << 328 K += s * (t+C); 337 } 329 } 338 K += C; 330 K += C; 339 } 331 } 340 else 332 else 341 { 333 { 342 G4double y = 1.0 / x / x; 334 G4double y = 1.0 / x / x; 343 K = 0.5 / x / pI0(x) * 335 K = 0.5 / x / pI0(x) * 344 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*( 336 (1.0 + y*(A0+y*(A1+y*(A2+y*(A3+y*(A4+y*(A5+y*(A6+y*A7)))))))); 345 } 337 } 346 338 347 return K; 339 return K; 348 } 340 } 349 ////////////////////////////////////////////// 341 //////////////////////////////////////////////////////////////////////////////// 350 // 342 // 351 G4double G4Bessel::pK1 (G4double x) 343 G4double G4Bessel::pK1 (G4double x) 352 { 344 { 353 G4double K = 0.0; 345 G4double K = 0.0; 354 if (x == 0.0) 346 if (x == 0.0) 355 K = 1.0E+307; 347 K = 1.0E+307; 356 else 348 else 357 K = (1.0/x - pI1(x)*pK0(x)) / pI0(x); 349 K = (1.0/x - pI1(x)*pK0(x)) / pI0(x); 358 return K; 350 return K; 359 } 351 } 360 ////////////////////////////////////////////// 352 //////////////////////////////////////////////////////////////////////////////// 361 // 353 // 362 354