Logo ROOT  
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
23#include <math.h>
24
25
26namespace ROOT {
27namespace Math {
28
29double KelvinFunctions::fgMin = 20;
30double KelvinFunctions::fgEpsilon = 1.e-20;
31
32double kSqrt2 = 1.4142135623730950488016887242097;
33double kPi = 3.14159265358979323846;
34double kEulerGamma = 0.577215664901532860606512090082402431042;
35
36
37/** \class KelvinFunctions
38This class calculates the Kelvin functions Ber(x), Bei(x), Ker(x),
39Kei(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
72double 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
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
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
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
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
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
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
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
487double 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
516double 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
547double 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
574double 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
603double 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
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
635double 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
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 double Phi(double x)
Utility function appearing in the asymptotic expansions of DKer(x) and DKei(x).
static double DBei(double x)
Calculates the first derivative of Bei(x).
static double M(double x)
Utility function appearing in the asymptotic expansions of DBer(x) and DBei(x).
static double DKer(double x)
Calculates the first derivative of Ker(x).
static double Ker(double x)
static double G1(double x)
Utility function appearing in the calculations of the Kelvin functions Bei(x) and Ber(x) (and their d...
static double DBer(double x)
Calculates the first derivative of Ber(x).
static double F2(double x)
Utility function appearing in the calculations of the Kelvin functions Kei(x) and Ker(x) (and their d...
static double Theta(double x)
Utility function appearing in the asymptotic expansions of DBer(x) and DBei(x).
static double N(double x)
Utility function appearing in the asymptotic expansions of DKer(x) and DKei(x).
static double Ber(double x)
static double F1(double x)
Utility function appearing in the calculations of the Kelvin functions Bei(x) and Ber(x) (and their d...
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 Bei(double x)
static double DKei(double x)
Calculates the first derivative of Kei(x).
static double Kei(double x)
double beta(double x, double y)
Calculates the beta function.
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double kEulerGamma
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static const double kSqrt2
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:96