Logo ROOT   6.14/05
Reference Guide
KelvinFunctions.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // The functions in this class have been imported by Jason Detwiler (jasondet@gmail.com) from
3 // CodeCogs GNU General Public License Agreement
4 // Copyright (C) 2004-2005 CodeCogs, Zyba Ltd, Broadwood, Holford, TA5 1DU,
5 // England.
6 //
7 // This program is free software; you can redistribute it and/or modify it
8 // under
9 // the terms of the GNU General Public License as published by CodeCogs.
10 // You must retain a copy of this licence in all copies.
11 //
12 // This program is distributed in the hope that it will be useful, but
13 // WITHOUT ANY
14 // WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A
16 // PARTICULAR PURPOSE. See the Adapted GNU General Public License for more
17 // details.
18 //
19 // *** THIS SOFTWARE CAN NOT BE USED FOR COMMERCIAL GAIN. ***
20 // ---------------------------------------------------------------------------------
21 
22 #include "Math/KelvinFunctions.h"
23 #include <math.h>
24 
25 
26 namespace ROOT {
27 namespace Math {
28 
29 double KelvinFunctions::fgMin = 20;
30 double KelvinFunctions::fgEpsilon = 1.e-20;
31 
32 double kSqrt2 = 1.4142135623730950488016887242097;
33 double kPi = 3.14159265358979323846;
34 double kEulerGamma = 0.577215664901532860606512090082402431042;
35 
36 
37 /** \class KelvinFunctions
38 This class calculates the Kelvin functions Ber(x), Bei(x), Ker(x),
39 Kei(x), and their first derivatives.
40 */
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 /// \f[
44 /// Ber(x) = Ber_{0}(x) = Re\left[J_{0}\left(x e^{3\pi i/4}\right)\right]
45 /// \f]
46 /// where x is real, and \f$J_{0}(z)\f$ is the zeroth-order Bessel
47 /// function of the first kind.
48 ///
49 /// If x < fgMin (=20), Ber(x) is computed according to its polynomial
50 /// approximation
51 /// \f[
52 /// Ber(x) = 1 + \sum_{n \geq 1}\frac{(-1)^{n}(x/2)^{4n}}{[(2n)!]^{2}}
53 /// \f]
54 /// For x > fgMin, Ber(x) is computed according to its asymptotic
55 /// expansion:
56 /// \f[
57 /// Ber(x) = \frac{e^{x/\sqrt{2}}}{\sqrt{2\pi x}} [F1(x) cos\alpha + G1(x) sin\alpha] - \frac{1}{\pi}Kei(x)
58 /// \f]
59 /// where \f$\alpha = \frac{x}{\sqrt{2}} - \frac{\pi}{8}\f$.
60 ///
61 /// See also F1() and G1().
62 ///
63 /// Begin_Macro
64 /// {
65 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
66 /// TF1 *fBer = new TF1("fBer","ROOT::Math::KelvinFunctions::Ber(x)",-10,10);
67 /// fBer->Draw();
68 /// return c;
69 /// }
70 /// End_Macro
71 
72 double KelvinFunctions::Ber(double x)
73 {
74  if (fabs(x) < fgEpsilon) return 1;
75 
76  if (fabs(x) < fgMin) {
77  double sum, factorial = 1, n = 1;
78  double term = 1, x_factor = x * x * x * x * 0.0625;
79 
80  sum = 1;
81 
82  do {
83  factorial = 4 * n * n * (2 * n - 1) * (2 * n - 1);
84  term *= (-1) / factorial * x_factor;
85  sum += term;
86  n += 1;
87  if (n > 1000) break;
88  } while (fabs(term) > fgEpsilon * sum);
89 
90  return sum;
91  } else {
92  double alpha = x / kSqrt2 - kPi / 8;
93  double value = F1(x) * cos(alpha) + G1(x) * sin(alpha);
94  value *= exp(x / kSqrt2) / sqrt(2 * kPi * x);
95  value -= Kei(x) / kPi;
96  return value;
97  }
98 }
99 
100 ////////////////////////////////////////////////////////////////////////////////
101 /// \f[
102 /// Bei(x) = Bei_{0}(x) = Im\left[J_{0}\left(x e^{3\pi i/4}\right)\right]
103 /// \f]
104 /// where x is real, and \f$J_{0}(z)\f$ is the zeroth-order Bessel
105 /// function of the first kind.
106 ///
107 /// If x < fgMin (=20), Bei(x) is computed according to its polynomial
108 /// approximation
109 /// \f[
110 /// Bei(x) = \sum_{n \geq 0}\frac{(-1)^{n}(x/2)^{4n+2}}{[(2n+1)!]^{2}}
111 /// \f]
112 /// For x > fgMin, Bei(x) is computed according to its asymptotic
113 /// expansion:
114 /// \f[
115 /// Bei(x) = \frac{e^{x/\sqrt{2}}}{\sqrt{2\pi x}} [F1(x) sin\alpha + G1(x) cos\alpha] - \frac{1}{\pi}Ker(x)
116 /// \f]
117 /// where \f$\alpha = \frac{x}{\sqrt{2}} - \frac{\pi}{8}\f$.
118 ///
119 /// See also F1() and G1().
120 ///
121 /// Begin_Macro
122 /// {
123 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
124 /// TF1 *fBei = new TF1("fBei","ROOT::Math::KelvinFunctions::Bei(x)",-10,10);
125 /// fBei->Draw();
126 /// return c;
127 /// }
128 /// End_Macro
129 
130 double KelvinFunctions::Bei(double x)
131 {
132 
133  if (fabs(x) < fgEpsilon) return 0;
134 
135  if (fabs(x) < fgMin) {
136  double sum, factorial = 1, n = 1;
137  double term = x * x * 0.25, x_factor = term * term;
138 
139  sum = term;
140 
141  do {
142  factorial = 4 * n * n * (2 * n + 1) * (2 * n + 1);
143  term *= (-1) / factorial * x_factor;
144  sum += term;
145  n += 1;
146  if (n > 1000) break;
147  } while (fabs(term) > fgEpsilon * sum);
148 
149  return sum;
150  } else {
151  double alpha = x / kSqrt2 - kPi / 8;
152  double value = F1(x) * sin(alpha) + G1(x) * cos(alpha);
153  value *= exp(x / kSqrt2) / sqrt(2 * kPi * x);
154  value += Ker(x) / kPi;
155  return value;
156  }
157 }
158 
159 
160 ////////////////////////////////////////////////////////////////////////////////
161 /// \f[
162 /// Ker(x) = Ker_{0}(x) = Re\left[K_{0}\left(x e^{3\pi i/4}\right)\right]
163 /// \f]
164 /// where x is real, and \f$K_{0}(z)\f$ is the zeroth-order modified
165 /// Bessel function of the second kind.
166 ///
167 /// If x < fgMin (=20), Ker(x) is computed according to its polynomial
168 /// approximation
169 /// \f[
170 /// Ker(x) = -\left(ln \frac{|x|}{2} + \gamma\right) Ber(x) + \left(\frac{\pi}{4} - \delta\right) Bei(x) + \sum_{n \geq 0} \frac{(-1)^{n}}{[(2n)!]^{2}} H_{2n} \left(\frac{x}{2}\right)^{4n}
171 /// \f]
172 /// where \f$\gamma = 0.577215664...\f$ is the Euler-Mascheroni constant,
173 /// \f$\delta = \pi\f$ for x < 0 and is otherwise zero, and
174 /// \f[
175 /// H_{n} = \sum_{k = 1}^{n} \frac{1}{k}
176 /// \f]
177 /// For x > fgMin, Ker(x) is computed according to its asymptotic
178 /// expansion:
179 /// \f[
180 /// Ker(x) = \sqrt{\frac{\pi}{2x}} e^{-x/\sqrt{2}} [F2(x) cos\beta + G2(x) sin\beta]
181 /// \f]
182 /// where \f$\beta = \frac{x}{\sqrt{2}} + \frac{\pi}{8}\f$.
183 ///
184 /// See also F2() and G2().
185 ///
186 /// Begin_Macro
187 /// {
188 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
189 /// TF1 *fKer = new TF1("fKer","ROOT::Math::KelvinFunctions::Ker(x)",-10,10);
190 /// fKer->Draw();
191 /// return c;
192 /// }
193 /// End_Macro
194 
195 double KelvinFunctions::Ker(double x)
196 {
197  if (fabs(x) < fgEpsilon) return 1E+100;
198 
199  if (fabs(x) < fgMin) {
200  double term = 1, x_factor = x * x * x * x * 0.0625;
201  double factorial = 1, harmonic = 0, n = 1, sum;
202  double delta = 0;
203  if(x < 0) delta = kPi;
204 
205  sum = - (log(fabs(x) * 0.5) + kEulerGamma) * Ber(x) + (kPi * 0.25 - delta) * Bei(x);
206 
207  do {
208  factorial = 4 * n * n * (2 * n - 1) * (2 * n - 1);
209  term *= (-1) / factorial * x_factor;
210  harmonic += 1 / (2 * n - 1 ) + 1 / (2 * n);
211  sum += term * harmonic;
212  n += 1;
213  if (n > 1000) break;
214  } while (fabs(term * harmonic) > fgEpsilon * sum);
215 
216  return sum;
217  } else {
218  double beta = x / kSqrt2 + kPi / 8;
219  double value = F2(x) * cos(beta) - G2(x) * sin(beta);
220  value *= sqrt(kPi / (2 * x)) * exp(- x / kSqrt2);
221  return value;
222  }
223 }
224 
225 
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// \f[
229 /// Kei(x) = Kei_{0}(x) = Im\left[K_{0}\left(x e^{3\pi i/4}\right)\right]
230 /// \f]
231 /// where x is real, and \f$K_{0}(z)\f$ is the zeroth-order modified
232 /// Bessel function of the second kind.
233 ///
234 /// If x < fgMin (=20), Kei(x) is computed according to its polynomial
235 /// approximation
236 /// \f[
237 /// Kei(x) = -\left(ln \frac{x}{2} + \gamma\right) Bei(x) - \left(\frac{\pi}{4} - \delta\right) Ber(x) + \sum_{n \geq 0} \frac{(-1)^{n}}{[(2n)!]^{2}} H_{2n} \left(\frac{x}{2}\right)^{4n+2}
238 /// \f]
239 /// where \f$\gamma = 0.577215664...\f$ is the Euler-Mascheroni constant,
240 /// \f$\delta = \pi\f$ for x < 0 and is otherwise zero, and
241 /// \f[
242 /// H_{n} = \sum_{k = 1}^{n} \frac{1}{k}
243 /// \f]
244 /// For x > fgMin, Kei(x) is computed according to its asymptotic
245 /// expansion:
246 /// \f[
247 /// Kei(x) = - \sqrt{\frac{\pi}{2x}} e^{-x/\sqrt{2}} [F2(x) sin\beta + G2(x) cos\beta]
248 /// \f]
249 /// where \f$\beta = \frac{x}{\sqrt{2}} + \frac{\pi}{8}\f$.
250 ///
251 /// See also F2() and G2().
252 ///
253 /// Begin_Macro
254 /// {
255 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
256 /// TF1 *fKei = new TF1("fKei","ROOT::Math::KelvinFunctions::Kei(x)",-10,10);
257 /// fKei->Draw();
258 /// return c;
259 /// }
260 /// End_Macro
261 
262 double KelvinFunctions::Kei(double x)
263 {
264  if (fabs(x) < fgEpsilon) return (-0.25 * kPi);
265 
266  if (fabs(x) < fgMin) {
267  double term = x * x * 0.25, x_factor = term * term;
268  double factorial = 1, harmonic = 1, n = 1, sum;
269  double delta = 0;
270  if(x < 0) delta = kPi;
271 
272  sum = term - (log(fabs(x) * 0.5) + kEulerGamma) * Bei(x) - (kPi * 0.25 - delta) * Ber(x);
273 
274  do {
275  factorial = 4 * n * n * (2 * n + 1) * (2 * n + 1);
276  term *= (-1) / factorial * x_factor;
277  harmonic += 1 / (2 * n) + 1 / (2 * n + 1);
278  sum += term * harmonic;
279  n += 1;
280  if (n > 1000) break;
281  } while (fabs(term * harmonic) > fgEpsilon * sum);
282 
283  return sum;
284  } else {
285  double beta = x / kSqrt2 + kPi / 8;
286  double value = - F2(x) * sin(beta) - G2(x) * cos(beta);
287  value *= sqrt(kPi / (2 * x)) * exp(- x / kSqrt2);
288  return value;
289  }
290 }
291 
292 
293 
294 ////////////////////////////////////////////////////////////////////////////////
295 /// Calculates the first derivative of Ber(x).
296 ///
297 /// If x < fgMin (=20), DBer(x) is computed according to the derivative of
298 /// the polynomial approximation of Ber(x). Otherwise it is computed
299 /// according to its asymptotic expansion
300 /// \f[
301 /// \frac{d}{dx} Ber(x) = M cos\left(\theta - \frac{\pi}{4}\right)
302 /// \f]
303 /// See also M() and Theta().
304 ///
305 /// Begin_Macro
306 /// {
307 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
308 /// TF1 *fDBer = new TF1("fDBer","ROOT::Math::KelvinFunctions::DBer(x)",-10,10);
309 /// fDBer->Draw();
310 /// return c;
311 /// }
312 /// End_Macro
313 
314 double KelvinFunctions::DBer(double x)
315 {
316  if (fabs(x) < fgEpsilon) return 0;
317 
318  if (fabs(x) < fgMin) {
319  double sum, factorial = 1, n = 1;
320  double term = - x * x * x * 0.0625, x_factor = - term * x;
321 
322  sum = term;
323 
324  do {
325  factorial = 4 * n * (n + 1) * (2 * n + 1) * (2 * n + 1);
326  term *= (-1) / factorial * x_factor;
327  sum += term;
328  n += 1;
329  if (n > 1000) break;
330  } while (fabs(term) > fgEpsilon * sum);
331 
332  return sum;
333  }
334  else return (M(x) * sin(Theta(x) - kPi / 4));
335 }
336 
337 
338 
339 ////////////////////////////////////////////////////////////////////////////////
340 /// Calculates the first derivative of Bei(x).
341 ///
342 /// If x < fgMin (=20), DBei(x) is computed according to the derivative of
343 /// the polynomial approximation of Bei(x). Otherwise it is computed
344 /// according to its asymptotic expansion
345 /// \f[
346 /// \frac{d}{dx} Bei(x) = M sin\left(\theta - \frac{\pi}{4}\right)
347 /// \f]
348 /// See also M() and Theta().
349 ///
350 /// Begin_Macro
351 /// {
352 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
353 /// TF1 *fDBei = new TF1("fDBei","ROOT::Math::KelvinFunctions::DBei(x)",-10,10);
354 /// fDBei->Draw();
355 /// return c;
356 /// }
357 /// End_Macro
358 
359 double KelvinFunctions::DBei(double x)
360 {
361  if (fabs(x) < fgEpsilon) return 0;
362 
363  if (fabs(x) < fgMin) {
364  double sum, factorial = 1, n = 1;
365  double term = x * 0.5, x_factor = x * x * x * x * 0.0625;
366 
367  sum = term;
368 
369  do {
370  factorial = 4 * n * n * (2 * n - 1) * (2 * n + 1);
371  term *= (-1) * x_factor / factorial;
372  sum += term;
373  n += 1;
374  if (n > 1000) break;
375  } while (fabs(term) > fgEpsilon * sum);
376 
377  return sum;
378  }
379  else return (M(x) * cos(Theta(x) - kPi / 4));
380 }
381 
382 
383 
384 ////////////////////////////////////////////////////////////////////////////////
385 /// Calculates the first derivative of Ker(x).
386 ///
387 /// If x < fgMin (=20), DKer(x) is computed according to the derivative of
388 /// the polynomial approximation of Ker(x). Otherwise it is computed
389 /// according to its asymptotic expansion
390 /// \f[
391 /// \frac{d}{dx} Ker(x) = N cos\left(\phi - \frac{\pi}{4}\right)
392 /// \f]
393 /// See also N() and Phi().
394 ///
395 /// Begin_Macro
396 /// {
397 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
398 /// TF1 *fDKer = new TF1("fDKer","ROOT::Math::KelvinFunctions::DKer(x)",-10,10);
399 /// fDKer->Draw();
400 /// return c;
401 /// }
402 /// End_Macro
403 
404 double KelvinFunctions::DKer(double x)
405 {
406  if (fabs(x) < fgEpsilon) return -1E+100;
407 
408  if (fabs(x) < fgMin) {
409  double term = - x * x * x * 0.0625, x_factor = - term * x;
410  double factorial = 1, harmonic = 1.5, n = 1, sum;
411  double delta = 0;
412  if(x < 0) delta = kPi;
413 
414  sum = 1.5 * term - Ber(x) / x - (log(fabs(x) * 0.5) + kEulerGamma) * DBer(x) + (0.25 * kPi - delta) * DBei(x);
415 
416  do {
417  factorial = 4 * n * (n + 1) * (2 * n + 1) * (2 * n + 1);
418  term *= (-1) / factorial * x_factor;
419  harmonic += 1 / (2 * n + 1 ) + 1 / (2 * n + 2);
420  sum += term * harmonic;
421  n += 1;
422  if (n > 1000) break;
423  } while (fabs(term * harmonic) > fgEpsilon * sum);
424 
425  return sum;
426  }
427  else return N(x) * sin(Phi(x) - kPi / 4);
428 }
429 
430 
431 
432 ////////////////////////////////////////////////////////////////////////////////
433 /// Calculates the first derivative of Kei(x).
434 ///
435 /// If x < fgMin (=20), DKei(x) is computed according to the derivative of
436 /// the polynomial approximation of Kei(x). Otherwise it is computed
437 /// according to its asymptotic expansion
438 /// \f[
439 /// \frac{d}{dx} Kei(x) = N sin\left(\phi - \frac{\pi}{4}\right)
440 /// \f]
441 /// See also N() and Phi().
442 ///
443 /// Begin_Macro
444 /// {
445 /// TCanvas *c = new TCanvas("c","c",0,0,500,300);
446 /// TF1 *fDKei = new TF1("fDKei","ROOT::Math::KelvinFunctions::DKei(x)",-10,10);
447 /// fDKei->Draw();
448 /// return c;
449 /// }
450 /// End_Macro
451 
452 double KelvinFunctions::DKei(double x)
453 {
454  if (fabs(x) < fgEpsilon) return 0;
455 
456  if (fabs(x) < fgMin) {
457  double term = 0.5 * x, x_factor = x * x * x * x * 0.0625;
458  double factorial = 1, harmonic = 1, n = 1, sum;
459  double delta = 0;
460  if(x < 0) delta = kPi;
461 
462  sum = term - Bei(x) / x - (log(fabs(x) * 0.5) + kEulerGamma) * DBei(x) - (kPi * 0.25 - delta) * DBer(x);
463 
464  do {
465  factorial = 4 * n * n * (2 * n - 1) * (2 * n + 1);
466  term *= (-1) / factorial * x_factor;
467  harmonic += 1 / (2 * n) + 1 / (2 * n + 1);
468  sum += term * harmonic;
469  n += 1;
470  if (n > 1000) break;
471  } while (fabs(term * harmonic) > fgEpsilon * sum);
472 
473  return sum;
474  }
475  else return N(x) * cos(Phi(x) - kPi / 4);
476 }
477 
478 
479 
480 ////////////////////////////////////////////////////////////////////////////////
481 /// Utility function appearing in the calculations of the Kelvin
482 /// functions Bei(x) and Ber(x) (and their derivatives). F1(x) is given by
483 /// \f[
484 /// F1(x) = 1 + \sum_{n \geq 1} \frac{\prod_{m=1}^{n}(2m - 1)^{2}}{n! (8x)^{n}} cos\left(\frac{n\pi}{4}\right)
485 /// \f]
486 
487 double KelvinFunctions::F1(double x)
488 {
489  double sum, term;
490  double prod = 1, x_factor = 8 * x, factorial = 1, n = 2;
491 
492  sum = kSqrt2 / (16 * x);
493 
494  do {
495  factorial *= n;
496  prod *= (2 * n - 1) * (2 * n - 1);
497  x_factor *= 8 * x;
498  term = prod / (factorial * x_factor) * cos(0.25 * n * kPi);
499  sum += term;
500  n += 1;
501  if (n > 1000) break;
502  } while (fabs(term) > fgEpsilon * sum);
503 
504  sum += 1;
505 
506  return sum;
507 }
508 
509 ////////////////////////////////////////////////////////////////////////////////
510 /// Utility function appearing in the calculations of the Kelvin
511 /// functions Kei(x) and Ker(x) (and their derivatives). F2(x) is given by
512 /// \f[
513 /// F2(x) = 1 + \sum_{n \geq 1} (-1)^{n} \frac{\prod_{m=1}^{n}(2m - 1)^{2}}{n! (8x)^{n}} cos\left(\frac{n\pi}{4}\right)
514 /// \f]
515 
516 double KelvinFunctions::F2(double x)
517 {
518  double sum, term;
519  double prod = 1, x_factor = 8 * x, factorial = 1, n = 2;
520 
521  sum = kSqrt2 / (16 * x);
522 
523  do {
524  factorial *= - n;
525  prod *= (2 * n - 1) * (2 * n - 1);
526  x_factor *= 8 * x;
527  term = (prod / (factorial * x_factor)) * cos(0.25 * n * kPi);
528  sum += term;
529  n += 1;
530  if (n > 1000) break;
531  } while (fabs(term) > fgEpsilon * sum);
532 
533  sum += 1;
534 
535  return sum;
536 }
537 
538 
539 
540 ////////////////////////////////////////////////////////////////////////////////
541 /// Utility function appearing in the calculations of the Kelvin
542 /// functions Bei(x) and Ber(x) (and their derivatives). G1(x) is given by
543 /// \f[
544 /// G1(x) = \sum_{n \geq 1} \frac{\prod_{m=1}^{n}(2m - 1)^{2}}{n! (8x)^{n}} sin\left(\frac{n\pi}{4}\right)
545 /// \f]
546 
547 double KelvinFunctions::G1(double x)
548 {
549  double sum, term;
550  double prod = 1, x_factor = 8 * x, factorial = 1, n = 2;
551 
552  sum = kSqrt2 / (16 * x);
553 
554  do {
555  factorial *= n;
556  prod *= (2 * n - 1) * (2 * n - 1);
557  x_factor *= 8 * x;
558  term = prod / (factorial * x_factor) * sin(0.25 * n * kPi);
559  sum += term;
560  n += 1;
561  if (n > 1000) break;
562  } while (fabs(term) > fgEpsilon * sum);
563 
564  return sum;
565 }
566 
567 ////////////////////////////////////////////////////////////////////////////////
568 /// Utility function appearing in the calculations of the Kelvin
569 /// functions Kei(x) and Ker(x) (and their derivatives). G2(x) is given by
570 /// \f[
571 /// G2(x) = \sum_{n \geq 1} (-1)^{n} \frac{\prod_{m=1}^{n}(2m - 1)^{2}}{n! (8x)^{n}} sin\left(\frac{n\pi}{4}\right)
572 /// \f]
573 
574 double KelvinFunctions::G2(double x)
575 {
576  double sum, term;
577  double prod = 1, x_factor = 8 * x, factorial = 1, n = 2;
578 
579  sum = kSqrt2 / (16 * x);
580 
581  do {
582  factorial *= - n;
583  prod *= (2 * n - 1) * (2 * n - 1);
584  x_factor *= 8 * x;
585  term = prod / (factorial * x_factor) * sin(0.25 * n * kPi);
586  sum += term;
587  n += 1;
588  if (n > 1000) break;
589  } while (fabs(term) > fgEpsilon * sum);
590 
591  return sum;
592 }
593 
594 
595 
596 ////////////////////////////////////////////////////////////////////////////////
597 /// Utility function appearing in the asymptotic expansions of DBer(x) and
598 /// DBei(x). M(x) is given by
599 /// \f[
600 /// M(x) = \frac{e^{x/\sqrt{2}}}{\sqrt{2\pi x}}\left(1 + \frac{1}{8\sqrt{2} x} + \frac{1}{256 x^{2}} - \frac{399}{6144\sqrt{2} x^{3}} + O\left(\frac{1}{x^{4}}\right)\right)
601 /// \f]
602 
603 double KelvinFunctions::M(double x)
604 {
605  double value = 1 + 1 / (8 * kSqrt2 * x) + 1 / (256 * x * x) - 399 / (6144 * kSqrt2 * x * x * x);
606  value *= exp(x / kSqrt2) / sqrt(2 * kPi * x);
607  return value;
608 }
609 
610 
611 
612 ////////////////////////////////////////////////////////////////////////////////
613 /// Utility function appearing in the asymptotic expansions of DBer(x) and
614 /// DBei(x). \f$\theta(x)\f$ is given by
615 /// \f[
616 /// \theta(x) = \frac{x}{\sqrt{2}} - \frac{\pi}{8} - \frac{1}{8\sqrt{2} x} - \frac{1}{16 x^{2}} - \frac{25}{384\sqrt{2} x^{3}} + O\left(\frac{1}{x^{5}}\right)
617 /// \f]
618 
619 double KelvinFunctions::Theta(double x)
620 {
621  double value = x / kSqrt2 - kPi / 8;
622  value -= 1 / (8 * kSqrt2 * x) + 1 / (16 * x * x) + 25 / (384 * kSqrt2 * x * x * x);
623  return value;
624 }
625 
626 
627 
628 ////////////////////////////////////////////////////////////////////////////////
629 /// Utility function appearing in the asymptotic expansions of DKer(x) and
630 /// DKei(x). N(x) is given by
631 /// \f[
632 /// N(x) = \sqrt{\frac{\pi}{2x}} e^{-x/\sqrt{2}} \left(1 - \frac{1}{8\sqrt{2} x} + \frac{1}{256 x^{2}} + \frac{399}{6144\sqrt{2} x^{3}} + O\left(\frac{1}{x^{4}}\right)\right)
633 /// \f]
634 
635 double KelvinFunctions::N(double x)
636 {
637  double value = 1 - 1 / (8 * kSqrt2 * x) + 1 / (256 * x * x) + 399 / (6144 * kSqrt2 * x * x * x);
638  value *= exp(- x / kSqrt2) * sqrt(kPi / (2 * x));
639  return value;
640 }
641 
642 
643 
644 ////////////////////////////////////////////////////////////////////////////////
645 /// Utility function appearing in the asymptotic expansions of DKer(x) and
646 /// DKei(x). \f$\phi(x)\f$ is given by
647 /// \f[
648 /// \phi(x) = - \frac{x}{\sqrt{2}} - \frac{\pi}{8} + \frac{1}{8\sqrt{2} x} - \frac{1}{16 x^{2}} + \frac{25}{384\sqrt{2} x^{3}} + O\left(\frac{1}{x^{5}}\right)
649 /// \f]
650 
651 double KelvinFunctions::Phi(double x)
652 {
653  double value = - x / kSqrt2 - kPi / 8;
654  value += 1 / (8 * kSqrt2 * x) - 1 / (16 * x * x) + 25 / (384 * kSqrt2 * x * x * x);
655  return value;
656 }
657 
658 
659 } // namespace Math
660 } // namespace ROOT
661 
662 
663 
static long int sum(long int i)
Definition: Factory.cxx:2258
static double Phi(double x)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static double Theta(double x)
static double F1(double x)
static double Ber(double x)
static const double kSqrt2
double kEulerGamma
double cos(double)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
double beta(double x, double y)
Calculates the beta function.
static double Bei(double x)
Double_t x[n]
Definition: legend1.C:17
static double G2(double x)
double sin(double)
static double M(double x)
static double G1(double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static double DKer(double x)
static double DKei(double x)
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
static double Kei(double x)
static double F2(double x)
static double N(double x)
static double Ker(double x)
static double DBei(double x)
Namespace for new Math classes and functions.
static double DBer(double x)
double exp(double)
const Int_t n
Definition: legend1.C:16
double log(double)