The Legendre Functions and Legendre Polynomials are described in Abramowitz & Stegun, Chapter 8.  These functions are declared in the header file gsl_sf_legendre.h.

Legendre Polynomials
--------------------

 -- Function: double gsl_sf_legendre_P1 (double X)
 -- Function: double gsl_sf_legendre_P2 (double X)
 -- Function: double gsl_sf_legendre_P3 (double X)
 -- Function: int gsl_sf_legendre_P1_e (double X, gsl_sf_result * RESULT)
 -- Function: int gsl_sf_legendre_P2_e (double X, gsl_sf_result * RESULT)
 -- Function: int gsl_sf_legendre_P3_e (double X, gsl_sf_result * RESULT)
     These functions evaluate the Legendre polynomials $P_l(x)$ P_l(x) using explicit representations for l=1, 2, 3.

 -- Function: double gsl_sf_legendre_Pl (int L, double X)
 -- Function: int gsl_sf_legendre_Pl_e (int L, double X, gsl_sf_result * RESULT)
     These functions evaluate the Legendre polynomial $P_l(x)$ P_l(x) for a specific value of L, X subject to $l \ge 0$ l >= 0, $|x| \le 1$ |x| <= 1

 -- Function: int gsl_sf_legendre_Pl_array (int LMAX, double X, double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_Pl_deriv_array (int LMAX, double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])

     These functions compute arrays of Legendre polynomials P_l(x) and derivatives dP_l(x)/dx, for l = 0, \dots, lmax, $|x| \le 1$ |x| <= 1

 -- Function: double gsl_sf_legendre_Q0 (double X)
 -- Function: int gsl_sf_legendre_Q0_e (double X, gsl_sf_result * RESULT)
     These routines compute the Legendre function Q_0(x) for x > -1, $x \ne 1$ x != 1.

 -- Function: double gsl_sf_legendre_Q1 (double X)
 -- Function: int gsl_sf_legendre_Q1_e (double X, gsl_sf_result * RESULT)
     These routines compute the Legendre function Q_1(x) for x > -1, $x \ne 1$ x != 1.

 -- Function: double gsl_sf_legendre_Ql (int L, double X)
 -- Function: int gsl_sf_legendre_Ql_e (int L, double X, gsl_sf_result * RESULT)
     These routines compute the Legendre function Q_l(x) for x > -1, $x \ne 1$ x != 1 and $l \ge 0$ l >= 0.

Associated Legendre Polynomials and Spherical Harmonics
-------------------------------------------------------

The following functions compute the associated Legendre polynomials P_l^m(x) which are solutions of the differential equation

     (1 - x^2) d^2 P_l^m(x) / dx^2 P_l^m(x) - 2x d/dx P_l^m(x) +
     ( l(l+1) - m^2 / (1 - x^2) ) P_l^m(x) = 0

where the degree l and order m satisfy 0 \le l and 0 \le m \le l. The functions P_l^m(x) grow combinatorially with l and can overflow for l larger than about 150. Alternatively, one may calculate normalized associated Legendre polynomials. There are a number of different normalization conventions, and these functions can be stably computed up to degree and order 2700. The following normalizations are provided:

Schmidt semi-normalization
     Schmidt semi-normalized associated Legendre polynomials are often used in the magnetics community and are defined as

          S_l^0(x) = P_l^0(x)
          S_l^m(x) = (-1)^m \sqrt((2(l-m)! / (l+m)!)) P_l^m(x), m > 0

     The factor of (-1)^m is called the Condon-Shortley phase factor and can be excluded if desired by setting the parameter csphase = 1 in the functions below.

Spherical Harmonic Normalization
     The associated Legendre polynomials suitable for calculating spherical harmonics are defined as

          Y_l^m(x) = (-1)^m \sqrt((2l + 1) * (l-m)! / (4 \pi) / (l+m)!) P_l^m(x)

     where again the phase factor (-1)^m can be included or excluded if desired.

Full Normalization
     The fully normalized associated Legendre polynomials are defined as

          N_l^m(x) = (-1)^m \sqrt((l + 1/2) * (l-m)! / (l+m)!) P_l^m(x)

     and have the property

          \int_(-1)^1 ( N_l^m(x) )^2 dx = 1

The normalized associated Legendre routines below use a recurrence relation which is stable up to a degree and order of about 2700. Beyond this, the computed functions could suffer from underflow leading to incorrect results. Routines are provided to compute first and second derivatives dP_l^m(x)/dx and d^2 P_l^m(x)/dx^2 as well as their alternate versions d P_l^m(\cos\theta)/d\theta and d^2 P_l^m(\cos\theta)/d\theta^2. While there is a simple scaling relationship between the two forms, the derivatives involving \theta are heavily used in spherical harmonic expansions and so these routines are also provided.

In the functions below, a parameter of type gsl_sf_legendre_t specifies the type of normalization to use. The possible values are

GSL_SF_LEGENDRE_NONE
     This specifies the computation of the unnormalized associated Legendre polynomials P_l^m(x).

GSL_SF_LEGENDRE_SCHMIDT
     This specifies the computation of the Schmidt semi-normalized associated Legendre polynomials S_l^m(x).

GSL_SF_LEGENDRE_SPHARM
     This specifies the computation of the spherical harmonic associated Legendre polynomials Y_l^m(x).

GSL_SF_LEGENDRE_FULL
     This specifies the computation of the fully normalized associated Legendre polynomials N_l^m(x).

 -- Function: int gsl_sf_legendre_array (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_array_e (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, const double CSPHASE, double RESULT_ARRAY[])
     These functions calculate all normalized associated Legendre polynomials for 0 \le l \le lmax and 0 \le m \le l for $|x| \le 1$ |x| <= 1. The NORM parameter specifies which normalization is used. The normalized P_l^m(x) values are stored in RESULT_ARRAY, whose minimum size can be obtained from calling gsl_sf_legendre_array_n. The array index of P_l^m(x) is obtained from calling gsl_sf_legendre_array_index(l, m). To include or exclude the Condon-Shortley phase factor of (-1)^m, set the parameter CSPHASE to either -1 or 1 respectively in the _e function. This factor is included by default.

 -- Function: int gsl_sf_legendre_deriv_array (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
 -- Function: int gsl_sf_legendre_deriv_array_e (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, const double CSPHASE, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     These functions calculate all normalized associated Legendre functions and their first derivatives up to degree LMAX for $|x| < 1$ |x| < 1. The parameter NORM specifies the normalization used. The normalized P_l^m(x) values and their derivatives dP_l^m(x)/dx are stored in RESULT_ARRAY and RESULT_DERIV_ARRAY respectively. To include or exclude the Condon-Shortley phase factor of (-1)^m, set the parameter CSPHASE to either -1 or 1 respectively in the _e function. This factor is included by default.

 -- Function: int gsl_sf_legendre_deriv_alt_array (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
 -- Function: int gsl_sf_legendre_deriv_alt_array_e (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, const double CSPHASE, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     These functions calculate all normalized associated Legendre functions and their (alternate) first derivatives up to degree LMAX for $|x| < 1$ |x| < 1. The normalized P_l^m(x) values and their derivatives dP_l^m(\cos\theta)/d\theta are stored in RESULT_ARRAY and RESULT_DERIV_ARRAY respectively. To include or exclude the Condon-Shortley phase factor of (-1)^m, set the parameter CSPHASE to either -1 or 1 respectively in the _e function. This factor is included by default.

 -- Function: int gsl_sf_legendre_deriv2_array (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[], double RESULT_DERIV2_ARRAY[])
 -- Function: int gsl_sf_legendre_deriv2_array_e (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, const double CSPHASE, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[], double RESULT_DERIV2_ARRAY[])
     These functions calculate all normalized associated Legendre functions and their first and second derivatives up to degree LMAX for $|x| < 1$ |x| < 1. The parameter NORM specifies the normalization used. The normalized P_l^m(x), their first derivatives dP_l^m(x)/dx, and their second derivatives d^2 P_l^m(x)/dx^2 are stored in RESULT_ARRAY, RESULT_DERIV_ARRAY, and RESULT_DERIV2_ARRAY respectively. To include or exclude the Condon-Shortley phase factor of (-1)^m, set the parameter CSPHASE to either -1 or 1 respectively in the _e function. This factor is included by default.

 -- Function: int gsl_sf_legendre_deriv2_alt_array (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[], double RESULT_DERIV2_ARRAY[])
 -- Function: int gsl_sf_legendre_deriv2_alt_array_e (const gsl_sf_legendre_t NORM, const size_t LMAX, const double X, const double CSPHASE, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[], double RESULT_DERIV2_ARRAY[])
     These functions calculate all normalized associated Legendre functions and their (alternate) first and second derivatives up to degree LMAX for $|x| < 1$ |x| < 1. The parameter NORM specifies the normalization used. The normalized P_l^m(x), their first derivatives dP_l^m(\cos\theta)/d\theta, and their second derivatives d^2 P_l^m(\cos\theta)/d\theta^2 are stored in RESULT_ARRAY, RESULT_DERIV_ARRAY, and RESULT_DERIV2_ARRAY respectively. To include or exclude the Condon-Shortley phase factor of (-1)^m, set the parameter CSPHASE to either -1 or 1 respectively in the _e function. This factor is included by default.

 -- Function: size_t gsl_sf_legendre_array_n (const size_t LMAX)
     This function returns the minimum array size for maximum degree LMAX needed for the array versions of the associated Legendre functions. Size is calculated as the total number of P_l^m(x) functions, plus extra space for precomputing multiplicative factors used in the recurrence relations.

 -- Function: size_t gsl_sf_legendre_array_index (const size_t L, const size_t M)
     This function returns the index into RESULT_ARRAY, RESULT_DERIV_ARRAY, or RESULT_DERIV2_ARRAY corresponding to P_l^m(x), P_l^'m(x), or P_l^''m(x). The index is given by l(l+1)/2 + m.

 -- Function: double gsl_sf_legendre_Plm (int L, int M, double X)
 -- Function: int gsl_sf_legendre_Plm_e (int L, int M, double X, gsl_sf_result * RESULT)
     These routines compute the associated Legendre polynomial P_l^m(x) for $m \ge 0$ m >= 0, $l \ge m$ l >= m, $|x| \le 1$ |x| <= 1.

 -- Function: double gsl_sf_legendre_sphPlm (int L, int M, double X)
 -- Function: int gsl_sf_legendre_sphPlm_e (int L, int M, double X, gsl_sf_result * RESULT)
     These routines compute the normalized associated Legendre polynomial $\sqrt(2l+1)/(4\pi) \sqrt(l-m)!/(l+m)! P_l^m(x)$ \sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x) suitable for use in spherical harmonics.  The parameters must satisfy $m \ge 0$ m >= 0, $l \ge m$ l >= m, $|x| \le 1$ |x| <= 1. Theses routines avoid the overflows that occur for the standard normalization of P_l^m(x).

 -- Function: int gsl_sf_legendre_Plm_array (int LMAX, int M, double X, double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_Plm_deriv_array (int LMAX, int M, double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     These functions are now deprecated and will be removed in a future release; see gsl_sf_legendre_array and gsl_sf_legendre_deriv_array.

 -- Function: int gsl_sf_legendre_sphPlm_array (int LMAX, int M, double X, double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_sphPlm_deriv_array (int LMAX, int M, double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     These functions are now deprecated and will be removed in a future release; see gsl_sf_legendre_array and gsl_sf_legendre_deriv_array.

 -- Function: int gsl_sf_legendre_array_size (const int LMAX, const int M)
     This function is now deprecated and will be removed in a future release.

Conical Functions
-----------------

The Conical Functions $P^\mu_-(1/2)+i\lambda(x)$ P^\mu_{-(1/2)+i\lambda}(x) and $Q^\mu_-(1/2)+i\lambda$ Q^\mu_{-(1/2)+i\lambda} are described in Abramowitz & Stegun, Section 8.12.

 -- Function: double gsl_sf_conicalP_half (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_half_e (double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the irregular Spherical Conical Function $P^1/2_-1/2 + i \lambda(x)$ P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_mhalf (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_mhalf_e (double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the regular Spherical Conical Function $P^-1/2_-1/2 + i \lambda(x)$ P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_0 (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_0_e (double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the conical function $P^0_-1/2 + i \lambda(x)$ P^0_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_1 (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_1_e (double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the conical function $P^1_-1/2 + i \lambda(x)$ P^1_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_sph_reg (int L, double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_sph_reg_e (int L, double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the Regular Spherical Conical Function $P^-1/2-l_-1/2 + i \lambda(x)$ P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, $l \ge -1$ l >= -1.

 -- Function: double gsl_sf_conicalP_cyl_reg (int M, double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_cyl_reg_e (int M, double LAMBDA, double X, gsl_sf_result * RESULT)
     These routines compute the Regular Cylindrical Conical Function $P^-m_-1/2 + i \lambda(x)$ P^{-m}_{-1/2 + i \lambda}(x) for x > -1, $m \ge -1$ m >= -1.

Radial Functions for Hyperbolic Space
-------------------------------------

The following spherical functions are specializations of Legendre functions which give the regular eigenfunctions of the Laplacian on a 3-dimensional hyperbolic space H3d.  Of particular interest is the flat limit, \lambda \to \infty, \eta \to 0, \lambda\eta fixed.

 -- Function: double gsl_sf_legendre_H3d_0 (double LAMBDA, double ETA)
 -- Function: int gsl_sf_legendre_H3d_0_e (double LAMBDA, double ETA, gsl_sf_result * RESULT)
     These routines compute the zeroth radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space, $$L^H3d_0(\lambda,\eta) := \sin(\lambda\eta) \over \lambda\sinh(\eta)$$ L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta)) for $\eta \ge 0$ \eta >= 0. In the flat limit this takes the form $L^H3d_0(\lambda,\eta) = j_0(\lambda\eta)$ L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta).

 -- Function: double gsl_sf_legendre_H3d_1 (double LAMBDA, double ETA)
 -- Function: int gsl_sf_legendre_H3d_1_e (double LAMBDA, double ETA, gsl_sf_result * RESULT)
     These routines compute the first radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space, $$L^H3d_1(\lambda,\eta) := 1\over\sqrt\lambda^2 + 1 \left(\sin(\lambda \eta)\over \lambda \sinh(\eta)\right) \left(\coth(\eta) - \lambda \cot(\lambda\eta)\right)$$ L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1} \sin(\lambda \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda \cot(\lambda\eta)) for $\eta \ge 0$ \eta >= 0. In the flat limit this takes the form $L^H3d_1(\lambda,\eta) = j_1(\lambda\eta)$ L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta).

 -- Function: double gsl_sf_legendre_H3d (int L, double LAMBDA, double ETA)
 -- Function: int gsl_sf_legendre_H3d_e (int L, double LAMBDA, double ETA, gsl_sf_result * RESULT)
     These routines compute the L-th radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space $\eta \ge 0$ \eta >= 0, $l \ge 0$ l >= 0. In the flat limit this takes the form $L^H3d_l(\lambda,\eta) = j_l(\lambda\eta)$ L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).

 -- Function: int gsl_sf_legendre_H3d_array (int LMAX, double LAMBDA, double ETA, double RESULT_ARRAY[])
     This function computes an array of radial eigenfunctions $L^H3d_l( \lambda, \eta)$ L^{H3d}_l(\lambda, \eta) for $0 \le l \le lmax$ 0 <= l <= lmax.