SpecFuncMathMore.cxx

Go to the documentation of this file.

// @(#)root/mathmore:$Id: SpecFuncMathMore.cxx 37160 2010-12-01 21:52:04Z moneta $
// Authors: L. Moneta, A. Zsenei   08/2005 

// Authors: Andras Zsenei & Lorenzo Moneta   06/2005 

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 , LCG ROOT MathLib Team                         *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

#include <cmath>

#ifndef PI
#define PI       3.14159265358979323846264338328      /* pi */
#endif


#include "gsl/gsl_sf_bessel.h"
#include "gsl/gsl_sf_legendre.h"
#include "gsl/gsl_sf_laguerre.h"
#include "gsl/gsl_sf_hyperg.h"
#include "gsl/gsl_sf_ellint.h"
//#include "gsl/gsl_sf_gamma.h"
#include "gsl/gsl_sf_expint.h"
#include "gsl/gsl_sf_zeta.h"
#include "gsl/gsl_sf_airy.h"
#include "gsl/gsl_sf_coupling.h"


namespace ROOT {
namespace Math {




// [5.2.1.1] associated Laguerre polynomials
// (26.x.12)

double assoc_laguerre(unsigned n, double m, double x) {
   
   return gsl_sf_laguerre_n(n, m, x);
   
}



// [5.2.1.2] associated Legendre functions
// (26.x.8)

double assoc_legendre(unsigned l, unsigned m, double x) {
   // add  (-1)^m
   return (m%2 == 0) ? gsl_sf_legendre_Plm(l, m, x) : -gsl_sf_legendre_Plm(l, m, x);
   
}





// [5.2.1.4] (complete) elliptic integral of the first kind
// (26.x.15.2) 

double comp_ellint_1(double k) {
   
   return gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
   
}



// [5.2.1.5] (complete) elliptic integral of the second kind
// (26.x.16.2)

double comp_ellint_2(double k) {
   
   return gsl_sf_ellint_Ecomp(k, GSL_PREC_DOUBLE);
   
}



// [5.2.1.6] (complete) elliptic integral of the third kind
// (26.x.17.2)
/**

There are two different definitions used for the elliptic
integral of the third kind:

P(\phi,k,n) = \int_0^\phi dt 1/((1 + nu \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

and

P(\phi,k,n) = \int_0^\phi dt 1/((1 - nu \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

the former is adopted by

- GSL 
     http://www.gnu.org/software/gsl/manual/gsl-ref_7.html#SEC95

- Planetmath
     http://planetmath.org/encyclopedia/EllipticIntegralsAndJacobiEllipticFunctions.html

- CERNLIB
     http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c346/top.html

   while the latter is used by

- Abramowitz and Stegun

- Mathematica
     http://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html

- C++ standard
     http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2004/n1687.pdf

   in order to be C++ compliant, we decided to use the latter, hence the 
   change of the sign in the function call to GSL.

   */

double comp_ellint_3(double n, double k) {
   
   return gsl_sf_ellint_P(PI/2.0, k, -n, GSL_PREC_DOUBLE);
   
}



// [5.2.1.7] confluent hypergeometric functions
// (26.x.14)

double conf_hyperg(double a, double b, double z) {
   
   return gsl_sf_hyperg_1F1(a, b, z);
   
}

//  confluent hypergeometric functions of second type

double conf_hypergU(double a, double b, double z) {
   
   return gsl_sf_hyperg_U(a, b, z);
   
}



// [5.2.1.8] regular modified cylindrical Bessel functions
// (26.x.4.1)

double cyl_bessel_i(double nu, double x) {
   
   return gsl_sf_bessel_Inu(nu, x);
   
}



// [5.2.1.9] cylindrical Bessel functions (of the first kind)
// (26.x.2)

double cyl_bessel_j(double nu, double x) {
   
   return gsl_sf_bessel_Jnu(nu, x);
   
}



// [5.2.1.10] irregular modified cylindrical Bessel functions
// (26.x.4.2)

double cyl_bessel_k(double nu, double x) {
   
   return gsl_sf_bessel_Knu(nu, x);
   
}



// [5.2.1.11] cylindrical Neumann functions;
// cylindrical Bessel functions (of the second kind)
// (26.x.3)

double cyl_neumann(double nu, double x) {
   
   return gsl_sf_bessel_Ynu(nu, x);
   
}



// [5.2.1.12] (incomplete) elliptic integral of the first kind
// phi in radians
// (26.x.15.1) 

double ellint_1(double k, double phi) {
   
   return gsl_sf_ellint_F(phi, k, GSL_PREC_DOUBLE);
   
}



// [5.2.1.13] (incomplete) elliptic integral of the second kind
// phi in radians
// (26.x.16.1)

double ellint_2(double k, double phi) {
   
   return gsl_sf_ellint_E(phi, k, GSL_PREC_DOUBLE);
   
}



// [5.2.1.14] (incomplete) elliptic integral of the third kind
// phi in radians
// (26.x.17.1)
/**

There are two different definitions used for the elliptic
integral of the third kind:

P(\phi,k,n) = \int_0^\phi dt 1/((1 + nu \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

and

P(\phi,k,n) = \int_0^\phi dt 1/((1 - nu \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

the former is adopted by

- GSL 
     http://www.gnu.org/software/gsl/manual/gsl-ref_7.html#SEC95

- Planetmath
     http://planetmath.org/encyclopedia/EllipticIntegralsAndJacobiEllipticFunctions.html

- CERNLIB
     http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c346/top.html

   while the latter is used by

- Abramowitz and Stegun

- Mathematica
     http://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html

- C++ standard
     http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2004/n1687.pdf

   in order to be C++ compliant, we decided to use the latter, hence the 
   change of the sign in the function call to GSL.

   */

double ellint_3(double n, double k, double phi) {
   
   return gsl_sf_ellint_P(phi, k, -n, GSL_PREC_DOUBLE);
   
}



// [5.2.1.15] exponential integral
// (26.x.20)

double expint(double x) {
   
   return gsl_sf_expint_Ei(x);
   
}



// [5.2.1.16] Hermite polynomials
// (26.x.10)

//double hermite(unsigned n, double x) {
//}




// [5.2.1.17] hypergeometric functions
// (26.x.13) 

double hyperg(double a, double b, double c, double x) {
   
   return gsl_sf_hyperg_2F1(a, b, c, x);
   
}



// [5.2.1.18] Laguerre polynomials
// (26.x.11)

double laguerre(unsigned n, double x) {
   return gsl_sf_laguerre_n(n, 0, x);
}




// [5.2.1.19] Legendre polynomials
// (26.x.7)

double legendre(unsigned l, double x) {
   
   return gsl_sf_legendre_Pl(l, x);
   
}



// [5.2.1.20] Riemann zeta function
// (26.x.22)

double riemann_zeta(double x) {
   
   return gsl_sf_zeta(x);
   
}



// [5.2.1.21] spherical Bessel functions of the first kind
// (26.x.5)

double sph_bessel(unsigned n, double x) {
   
   return gsl_sf_bessel_jl(n, x);
   
}



// [5.2.1.22] spherical associated Legendre functions
// (26.x.9) ??????????

double sph_legendre(unsigned l, unsigned m, double theta) {

   return gsl_sf_legendre_sphPlm(l, m, std::cos(theta));

}




// [5.2.1.23] spherical Neumann functions
// (26.x.6)

double sph_neumann(unsigned n, double x) {
   
   return gsl_sf_bessel_yl(n, x);
   
} 

// Airy function Ai

double airy_Ai(double x) {
   
   return gsl_sf_airy_Ai(x, GSL_PREC_DOUBLE);
   
} 

// Airy function Bi

double airy_Bi(double x) {
   
   return gsl_sf_airy_Bi(x, GSL_PREC_DOUBLE);
   
}

// Derivative of the Airy function Ai

double airy_Ai_deriv(double x) {
   
   return gsl_sf_airy_Ai_deriv(x, GSL_PREC_DOUBLE);
   
}

// Derivative of the Airy function Bi

double airy_Bi_deriv(double x) {
   
   return gsl_sf_airy_Bi_deriv(x, GSL_PREC_DOUBLE);
   
}

// s-th zero of the Airy function Ai

double airy_zero_Ai(unsigned int s) {
   
   return gsl_sf_airy_zero_Ai(s);
   
}

// s-th zero of the Airy function Bi

double airy_zero_Bi(unsigned int s) {
   
   return gsl_sf_airy_zero_Bi(s);
   
}

// s-th zero of the derivative of the Airy function Ai

double airy_zero_Ai_deriv(unsigned int s) {
   
   return gsl_sf_airy_zero_Ai_deriv(s);
   
}

// s-th zero of the derivative of the Airy function Bi

double airy_zero_Bi_deriv(unsigned int s) {
   
   return gsl_sf_airy_zero_Bi_deriv(s);
   
}

// wigner coefficient 

double wigner_3j(int ja, int jb, int jc, int ma, int mb, int mc) {
   return gsl_sf_coupling_3j(ja,jb,jc,ma,mb,mc);
}

double wigner_6j(int ja, int jb, int jc, int jd, int je, int jf) { 
   return gsl_sf_coupling_6j(ja,jb,jc,jd,je,jf);
}

double wigner_9j(int ja, int jb, int jc, int jd, int je, int jf, int jg, int jh, int ji) { 
   return gsl_sf_coupling_9j(ja,jb,jc,jd,je,jf,jg,jh,ji);
}

} // namespace Math
} // namespace ROOT

---