Logo ROOT  
Reference Guide
SpecFuncMathMore.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4 // Authors: Andras Zsenei & Lorenzo Moneta 06/2005
5 
6 /**********************************************************************
7  * *
8  * Copyright (c) 2005 , LCG ROOT MathLib Team *
9  * *
10  * *
11  **********************************************************************/
12 
13 #include <cmath>
14 
15 #ifndef PI
16 #define PI 3.14159265358979323846264338328 /* pi */
17 #endif
18 
19 
20 #include "gsl/gsl_sf_bessel.h"
21 #include "gsl/gsl_sf_legendre.h"
22 #include "gsl/gsl_sf_laguerre.h"
23 #include "gsl/gsl_sf_hyperg.h"
24 #include "gsl/gsl_sf_ellint.h"
25 //#include "gsl/gsl_sf_gamma.h"
26 #include "gsl/gsl_sf_expint.h"
27 #include "gsl/gsl_sf_zeta.h"
28 #include "gsl/gsl_sf_airy.h"
29 #include "gsl/gsl_sf_coupling.h"
30 
31 
32 namespace ROOT {
33 namespace Math {
34 
35 
36 
37 
38 // [5.2.1.1] associated Laguerre polynomials
39 // (26.x.12)
40 
41 double assoc_laguerre(unsigned n, double m, double x) {
42 
43  return gsl_sf_laguerre_n(n, m, x);
44 
45 }
46 
47 
48 
49 // [5.2.1.2] associated Legendre functions
50 // (26.x.8)
51 
52 double assoc_legendre(unsigned l, unsigned m, double x) {
53  // add (-1)^m
54  return (m%2 == 0) ? gsl_sf_legendre_Plm(l, m, x) : -gsl_sf_legendre_Plm(l, m, x);
55 
56 }
57 
58 // Shortcut for RooFit to call the gsl legendre functions without the branches in the above implementation.
59 namespace internal{
60  double legendre(unsigned l, unsigned m, double x) {
61  return gsl_sf_legendre_Plm(l, m, x);
62  }
63 }
64 
65 
66 // [5.2.1.4] (complete) elliptic integral of the first kind
67 // (26.x.15.2)
68 
69 double comp_ellint_1(double k) {
70 
71  return gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
72 
73 }
74 
75 
76 
77 // [5.2.1.5] (complete) elliptic integral of the second kind
78 // (26.x.16.2)
79 
80 double comp_ellint_2(double k) {
81 
82  return gsl_sf_ellint_Ecomp(k, GSL_PREC_DOUBLE);
83 
84 }
85 
86 
87 
88 // [5.2.1.6] (complete) elliptic integral of the third kind
89 // (26.x.17.2)
90 /**
91 Complete elliptic integral of the third kind.
92 
93 There are two different definitions used for the elliptic
94 integral of the third kind:
95 
96 \f[
97 P(\phi,k,n) = \int_0^\phi \frac{dt}{(1 + n \sin^2{t})\sqrt{1 - k^2 \sin^2{t}}}
98 \f]
99 
100 and
101 
102 \f[
103 P(\phi,k,n) = \int_0^\phi \frac{dt}{(1 - n \sin^2{t})\sqrt{1 - k^2 \sin^2{t}}}
104 \f]
105 
106 the former is adopted by
107 
108 - GSL
109  http://www.gnu.org/software/gsl/manual/gsl-ref_7.html#SEC95
110 
111 - Planetmath
112  http://planetmath.org/encyclopedia/EllipticIntegralsAndJacobiEllipticFunctions.html
113 
114 - CERNLIB
115  https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/c346/top.html
116 
117  while the latter is used by
118 
119 - Abramowitz and Stegun
120 
121 - Mathematica
122  http://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html
123 
124 - C++ standard
125  http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2004/n1687.pdf
126 
127  in order to be C++ compliant, we decided to use the latter, hence the
128  change of the sign in the function call to GSL.
129 
130  */
131 
132 double comp_ellint_3(double n, double k) {
133 
134  return gsl_sf_ellint_P(PI/2.0, k, -n, GSL_PREC_DOUBLE);
135 
136 }
137 
138 
139 
140 // [5.2.1.7] confluent hypergeometric functions
141 // (26.x.14)
142 
143 double conf_hyperg(double a, double b, double z) {
144 
145  return gsl_sf_hyperg_1F1(a, b, z);
146 
147 }
148 
149 // confluent hypergeometric functions of second type
150 
151 double conf_hypergU(double a, double b, double z) {
152 
153  return gsl_sf_hyperg_U(a, b, z);
154 
155 }
156 
157 
158 
159 // [5.2.1.8] regular modified cylindrical Bessel functions
160 // (26.x.4.1)
161 
162 double cyl_bessel_i(double nu, double x) {
163 
164  return gsl_sf_bessel_Inu(nu, x);
165 
166 }
167 
168 
169 
170 // [5.2.1.9] cylindrical Bessel functions (of the first kind)
171 // (26.x.2)
172 
173 double cyl_bessel_j(double nu, double x) {
174 
175  return gsl_sf_bessel_Jnu(nu, x);
176 
177 }
178 
179 
180 
181 // [5.2.1.10] irregular modified cylindrical Bessel functions
182 // (26.x.4.2)
183 
184 double cyl_bessel_k(double nu, double x) {
185 
186  return gsl_sf_bessel_Knu(nu, x);
187 
188 }
189 
190 
191 
192 // [5.2.1.11] cylindrical Neumann functions;
193 // cylindrical Bessel functions (of the second kind)
194 // (26.x.3)
195 
196 double cyl_neumann(double nu, double x) {
197 
198  return gsl_sf_bessel_Ynu(nu, x);
199 
200 }
201 
202 
203 
204 // [5.2.1.12] (incomplete) elliptic integral of the first kind
205 // phi in radians
206 // (26.x.15.1)
207 
208 double ellint_1(double k, double phi) {
209 
210  return gsl_sf_ellint_F(phi, k, GSL_PREC_DOUBLE);
211 
212 }
213 
214 
215 
216 // [5.2.1.13] (incomplete) elliptic integral of the second kind
217 // phi in radians
218 // (26.x.16.1)
219 
220 double ellint_2(double k, double phi) {
221 
222  return gsl_sf_ellint_E(phi, k, GSL_PREC_DOUBLE);
223 
224 }
225 
226 
227 
228 // [5.2.1.14] (incomplete) elliptic integral of the third kind
229 // phi in radians
230 // (26.x.17.1)
231 /**
232 
233 Incomplete elliptic integral of the third kind.
234 
235 There are two different definitions used for the elliptic
236 integral of the third kind:
237 
238 \f[
239 P(\phi,k,n) = \int_0^\phi \frac{dt}{(1 + n \sin^2{t})\sqrt{1 - k^2 \sin^2{t}}}
240 \f]
241 
242 and
243 
244 \f[
245 P(\phi,k,n) = \int_0^\phi \frac{dt}{(1 - n \sin^2{t})\sqrt{1 - k^2 \sin^2{t}}}
246 \f]
247 
248 the former is adopted by
249 
250 - GSL
251  http://www.gnu.org/software/gsl/manual/gsl-ref_7.html#SEC95
252 
253 - Planetmath
254  http://planetmath.org/encyclopedia/EllipticIntegralsAndJacobiEllipticFunctions.html
255 
256 - CERNLIB
257  https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/c346/top.html
258 
259  while the latter is used by
260 
261 - Abramowitz and Stegun
262 
263 - Mathematica
264  http://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html
265 
266 - C++ standard
267  http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2004/n1687.pdf
268 
269  in order to be C++ compliant, we decided to use the latter, hence the
270  change of the sign in the function call to GSL.
271 
272  */
273 
274 double ellint_3(double n, double k, double phi) {
275 
276  return gsl_sf_ellint_P(phi, k, -n, GSL_PREC_DOUBLE);
277 
278 }
279 
280 
281 
282 // [5.2.1.15] exponential integral
283 // (26.x.20)
284 
285 double expint(double x) {
286 
287  return gsl_sf_expint_Ei(x);
288 
289 }
290 
291 
292 // Generalization of expint(x)
293 //
294 double expint_n(int n, double x) {
295 
296  return gsl_sf_expint_En(n, x);
297 
298 }
299 
300 
301 
302 // [5.2.1.16] Hermite polynomials
303 // (26.x.10)
304 
305 //double hermite(unsigned n, double x) {
306 //}
307 
308 
309 
310 
311 // [5.2.1.17] hypergeometric functions
312 // (26.x.13)
313 
314 double hyperg(double a, double b, double c, double x) {
315 
316  return gsl_sf_hyperg_2F1(a, b, c, x);
317 
318 }
319 
320 
321 
322 // [5.2.1.18] Laguerre polynomials
323 // (26.x.11)
324 
325 double laguerre(unsigned n, double x) {
326  return gsl_sf_laguerre_n(n, 0, x);
327 }
328 
329 
330 
331 
332 // [5.2.1.19] Legendre polynomials
333 // (26.x.7)
334 
335 double legendre(unsigned l, double x) {
336 
337  return gsl_sf_legendre_Pl(l, x);
338 
339 }
340 
341 
342 
343 // [5.2.1.20] Riemann zeta function
344 // (26.x.22)
345 
346 double riemann_zeta(double x) {
347 
348  return gsl_sf_zeta(x);
349 
350 }
351 
352 
353 
354 // [5.2.1.21] spherical Bessel functions of the first kind
355 // (26.x.5)
356 
357 double sph_bessel(unsigned n, double x) {
358 
359  return gsl_sf_bessel_jl(n, x);
360 
361 }
362 
363 
364 
365 // [5.2.1.22] spherical associated Legendre functions
366 // (26.x.9) ??????????
367 
368 double sph_legendre(unsigned l, unsigned m, double theta) {
369 
370  return gsl_sf_legendre_sphPlm(l, m, std::cos(theta));
371 
372 }
373 
374 
375 
376 
377 // [5.2.1.23] spherical Neumann functions
378 // (26.x.6)
379 
380 double sph_neumann(unsigned n, double x) {
381 
382  return gsl_sf_bessel_yl(n, x);
383 
384 }
385 
386 // Airy function Ai
387 
388 double airy_Ai(double x) {
389 
390  return gsl_sf_airy_Ai(x, GSL_PREC_DOUBLE);
391 
392 }
393 
394 // Airy function Bi
395 
396 double airy_Bi(double x) {
397 
398  return gsl_sf_airy_Bi(x, GSL_PREC_DOUBLE);
399 
400 }
401 
402 // Derivative of the Airy function Ai
403 
404 double airy_Ai_deriv(double x) {
405 
406  return gsl_sf_airy_Ai_deriv(x, GSL_PREC_DOUBLE);
407 
408 }
409 
410 // Derivative of the Airy function Bi
411 
412 double airy_Bi_deriv(double x) {
413 
414  return gsl_sf_airy_Bi_deriv(x, GSL_PREC_DOUBLE);
415 
416 }
417 
418 // s-th zero of the Airy function Ai
419 
420 double airy_zero_Ai(unsigned int s) {
421 
422  return gsl_sf_airy_zero_Ai(s);
423 
424 }
425 
426 // s-th zero of the Airy function Bi
427 
428 double airy_zero_Bi(unsigned int s) {
429 
430  return gsl_sf_airy_zero_Bi(s);
431 
432 }
433 
434 // s-th zero of the derivative of the Airy function Ai
435 
436 double airy_zero_Ai_deriv(unsigned int s) {
437 
438  return gsl_sf_airy_zero_Ai_deriv(s);
439 
440 }
441 
442 // s-th zero of the derivative of the Airy function Bi
443 
444 double airy_zero_Bi_deriv(unsigned int s) {
445 
446  return gsl_sf_airy_zero_Bi_deriv(s);
447 
448 }
449 
450 // wigner coefficient
451 
452 double wigner_3j(int ja, int jb, int jc, int ma, int mb, int mc) {
453  return gsl_sf_coupling_3j(ja,jb,jc,ma,mb,mc);
454 }
455 
456 double wigner_6j(int ja, int jb, int jc, int jd, int je, int jf) {
457  return gsl_sf_coupling_6j(ja,jb,jc,jd,je,jf);
458 }
459 
460 double wigner_9j(int ja, int jb, int jc, int jd, int je, int jf, int jg, int jh, int ji) {
461  return gsl_sf_coupling_9j(ja,jb,jc,jd,je,jf,jg,jh,ji);
462 }
463 
464 } // namespace Math
465 } // namespace ROOT
c
#define c(i)
Definition: RSha256.hxx:119
l
auto * l
Definition: textangle.C:4
m
auto * m
Definition: textangle.C:8
n
const Int_t n
Definition: legend1.C:16
PI
#define PI
Definition: SpecFuncMathMore.cxx:16
ROOT::Math::sph_legendre
double sph_legendre(unsigned l, unsigned m, double theta)
Computes the spherical (normalized) associated Legendre polynomials, or spherical harmonic without az...
Definition: SpecFuncMathMore.cxx:368
ROOT::Math::cyl_bessel_k
double cyl_bessel_k(double nu, double x)
Calculates the modified Bessel functions of the second kind (also called irregular modified (cylindri...
Definition: SpecFuncMathMore.cxx:184
ROOT::Math::cyl_bessel_j
double cyl_bessel_j(double nu, double x)
Calculates the (cylindrical) Bessel functions of the first kind (also called regular (cylindrical) Be...
Definition: SpecFuncMathMore.cxx:173
ROOT::Math::airy_zero_Bi
double airy_zero_Bi(unsigned int s)
Calculates the zeroes of the Airy function Bi.
Definition: SpecFuncMathMore.cxx:428
ROOT::Math::airy_zero_Ai_deriv
double airy_zero_Ai_deriv(unsigned int s)
Calculates the zeroes of the derivative of the Airy function Ai.
Definition: SpecFuncMathMore.cxx:436
ROOT::Math::ellint_2
double ellint_2(double k, double phi)
Calculates the complete elliptic integral of the second kind.
Definition: SpecFuncMathMore.cxx:220
ROOT::Math::wigner_6j
double wigner_6j(int two_ja, int two_jb, int two_jc, int two_jd, int two_je, int two_jf)
Calculates the Wigner 6j coupling coefficients.
Definition: SpecFuncMathMore.cxx:456
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:168
ROOT::Math::legendre
double legendre(unsigned l, double x)
Calculates the Legendre polynomials.
Definition: SpecFuncMathMore.cxx:335
x
Double_t x[n]
Definition: legend1.C:17
cos
double cos(double)
ROOT::Math::conf_hyperg
double conf_hyperg(double a, double b, double z)
Calculates the confluent hypergeometric functions of the first kind.
Definition: SpecFuncMathMore.cxx:143
b
#define b(i)
Definition: RSha256.hxx:118
ROOT::Math::comp_ellint_1
double comp_ellint_1(double k)
Calculates the complete elliptic integral of the first kind.
Definition: SpecFuncMathMore.cxx:69
ROOT::Math::internal::legendre
double legendre(unsigned l, unsigned m, double x)
Definition: SpecFuncMathMore.cxx:60
ROOT::Math::sph_bessel
double sph_bessel(unsigned n, double x)
Calculates the spherical Bessel functions of the first kind (also called regular spherical Bessel fun...
Definition: SpecFuncMathMore.cxx:357
ROOT::Math::airy_zero_Bi_deriv
double airy_zero_Bi_deriv(unsigned int s)
Calculates the zeroes of the derivative of the Airy function Bi.
Definition: SpecFuncMathMore.cxx:444
ROOT::Math::hyperg
double hyperg(double a, double b, double c, double x)
Calculates Gauss' hypergeometric function.
Definition: SpecFuncMathMore.cxx:314
ROOT::Math::comp_ellint_2
double comp_ellint_2(double k)
Calculates the complete elliptic integral of the second kind.
Definition: SpecFuncMathMore.cxx:80
ROOT::Math::wigner_9j
double wigner_9j(int two_ja, int two_jb, int two_jc, int two_jd, int two_je, int two_jf, int two_jg, int two_jh, int two_ji)
Calculates the Wigner 9j coupling coefficients.
Definition: SpecFuncMathMore.cxx:460
ROOT::Math::expint
double expint(double x)
Calculates the exponential integral.
Definition: SpecFuncMathMore.cxx:285
a
auto * a
Definition: textangle.C:12
ROOT::Math::assoc_legendre
double assoc_legendre(unsigned l, unsigned m, double x)
Computes the associated Legendre polynomials.
Definition: SpecFuncMathMore.cxx:52
ROOT::Math::ellint_3
double ellint_3(double n, double k, double phi)
Calculates the complete elliptic integral of the third kind.
Definition: SpecFuncMathMore.cxx:274
ROOT::Math::riemann_zeta
double riemann_zeta(double x)
Calculates the Riemann zeta function.
Definition: SpecFuncMathMore.cxx:346
ROOT::Math::airy_Bi_deriv
double airy_Bi_deriv(double x)
Calculates the derivative of the Airy function Bi.
Definition: SpecFuncMathMore.cxx:412
ROOT::Math::assoc_laguerre
double assoc_laguerre(unsigned n, double m, double x)
Computes the generalized Laguerre polynomials for .
Definition: SpecFuncMathMore.cxx:41
ROOT::Math::airy_Ai_deriv
double airy_Ai_deriv(double x)
Calculates the derivative of the Airy function Ai.
Definition: SpecFuncMathMore.cxx:404
ROOT::Math::ellint_1
double ellint_1(double k, double phi)
Calculates the incomplete elliptic integral of the first kind.
Definition: SpecFuncMathMore.cxx:208
ROOT::Math::airy_zero_Ai
double airy_zero_Ai(unsigned int s)
Calculates the zeroes of the Airy function Ai.
Definition: SpecFuncMathMore.cxx:420
ROOT::Math::expint_n
double expint_n(int n, double x)
Definition: SpecFuncMathMore.cxx:294
ROOT::Math::airy_Bi
double airy_Bi(double x)
Calculates the Airy function Bi.
Definition: SpecFuncMathMore.cxx:396
ROOT::Math::laguerre
double laguerre(unsigned n, double x)
Calculates the Laguerre polynomials.
Definition: SpecFuncMathMore.cxx:325
ROOT::Math::comp_ellint_3
double comp_ellint_3(double n, double k)
Calculates the complete elliptic integral of the third kind.
Definition: SpecFuncMathMore.cxx:132
ROOT::Math::wigner_3j
double wigner_3j(int two_ja, int two_jb, int two_jc, int two_ma, int two_mb, int two_mc)
Calculates the Wigner 3j coupling coefficients.
Definition: SpecFuncMathMore.cxx:452
ROOT::Math::sph_neumann
double sph_neumann(unsigned n, double x)
Calculates the spherical Bessel functions of the second kind (also called irregular spherical Bessel ...
Definition: SpecFuncMathMore.cxx:380
ROOT::Math::cyl_neumann
double cyl_neumann(double nu, double x)
Calculates the (cylindrical) Bessel functions of the second kind (also called irregular (cylindrical)...
Definition: SpecFuncMathMore.cxx:196
ROOT::Math::cyl_bessel_i
double cyl_bessel_i(double nu, double x)
Calculates the modified Bessel function of the first kind (also called regular modified (cylindrical)...
Definition: SpecFuncMathMore.cxx:162
ROOT::Math::airy_Ai
double airy_Ai(double x)
Calculates the Airy function Ai.
Definition: SpecFuncMathMore.cxx:388
ROOT::Math::conf_hypergU
double conf_hypergU(double a, double b, double z)
Calculates the confluent hypergeometric functions of the second kind, known also as Kummer function o...
Definition: SpecFuncMathMore.cxx:151
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Math