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