ROOT   6.10/09 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
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 ///
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 ///
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 ///
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 ///
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]
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]
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]
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]
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:2162
static double Ber(double x)
where x is real, and is the zeroth-order Bessel function of the first kind.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static double DBei(double x)
Calculates the first derivative of Bei(x).
static const double kSqrt2
double kEulerGamma
Float_t delta
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 F1(double x)
Utility function appearing in the calculations of the Kelvin functions Bei(x) and Ber(x) (and their d...
Double_t x[n]
Definition: legend1.C:17
static double G1(double x)
Utility function appearing in the calculations of the Kelvin functions Bei(x) and Ber(x) (and their d...
double sin(double)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static double Kei(double x)
where x is real, and is the zeroth-order modified Bessel function of the second kind...
static double DKer(double x)
Calculates the first derivative of Ker(x).
static double M(double x)
Utility function appearing in the asymptotic expansions of DBer(x) and DBei(x).
constexpr Double_t E()
Definition: TMath.h:74
static double Bei(double x)
where x is real, and is the zeroth-order Bessel function of the first kind.
static double Theta(double x)
Utility function appearing in the asymptotic expansions of DBer(x) and DBei(x).
Namespace for new Math classes and functions.
static double G2(double x)
Utility function appearing in the calculations of the Kelvin functions Kei(x) and Ker(x) (and their d...
static double DBer(double x)
Calculates the first derivative of Ber(x).
static double Phi(double x)
Utility function appearing in the asymptotic expansions of DKer(x) and DKei(x).
static double Ker(double x)
where x is real, and is the zeroth-order modified Bessel function of the second kind...
static double DKei(double x)
Calculates the first derivative of Kei(x).
double exp(double)
static double F2(double x)
Utility function appearing in the calculations of the Kelvin functions Kei(x) and Ker(x) (and their d...
const Int_t n
Definition: legend1.C:16
double log(double)
static double N(double x)
Utility function appearing in the asymptotic expansions of DKer(x) and DKei(x).