The routines described in this section compute the angular and radial Mathieu functions, and their characteristic values.  Mathieu functions are the solutions of the following two differential equations:

     d^2y/dv^2 + (a - 2q\cos 2v)y = 0
     d^2f/du^2 - (a - 2q\cosh 2u)f = 0

The angular Mathieu functions ce_r(x,q), se_r(x,q) are the even and odd periodic solutions of the first equation, which is known as Mathieu's equation. These exist only for the discrete sequence of  characteristic values a=a_r(q) (even-periodic) and a=b_r(q) (odd-periodic).

The radial Mathieu functions $Mc^(j)_r(z,q)$ Mc^{(j)}_{r}(z,q), $Ms^(j)_r(z,q)$ Ms^{(j)}_{r}(z,q) are the solutions of the second equation, which is referred to as Mathieu's modified equation.  The radial Mathieu functions of the first, second, third and fourth kind are denoted by the parameter j, which takes the value 1, 2, 3 or 4.

For more information on the Mathieu functions, see Abramowitz and Stegun, Chapter 20.  These functions are defined in the header file gsl_sf_mathieu.h.

Mathieu Function Workspace
--------------------------

The Mathieu functions can be computed for a single order or for multiple orders, using array-based routines.  The array-based routines require a preallocated workspace.

 -- Function: gsl_sf_mathieu_workspace * gsl_sf_mathieu_alloc (size_t N, double QMAX)
     This function returns a workspace for the array versions of the Mathieu routines.  The arguments N and QMAX specify the maximum order and q-value of Mathieu functions which can be computed with this workspace.

 -- Function: void gsl_sf_mathieu_free (gsl_sf_mathieu_workspace * WORK)
     This function frees the workspace WORK.

Mathieu Function Characteristic Values
--------------------------------------

 -- Function: int gsl_sf_mathieu_a (int N, double Q)
 -- Function: int gsl_sf_mathieu_a_e (int N, double Q, gsl_sf_result * RESULT)
 -- Function: int gsl_sf_mathieu_b (int N, double Q)
 -- Function: int gsl_sf_mathieu_b_e (int N, double Q, gsl_sf_result * RESULT)
     These routines compute the characteristic values a_n(q), b_n(q) of the Mathieu functions ce_n(q,x) and se_n(q,x), respectively.

 -- Function: int gsl_sf_mathieu_a_array (int ORDER_MIN, int ORDER_MAX, double Q, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_b_array (int ORDER_MIN, int ORDER_MAX, double Q, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
     These routines compute a series of Mathieu characteristic values a_n(q), b_n(q) for n from ORDER_MIN to ORDER_MAX inclusive, storing the results in the array RESULT_ARRAY.

Angular Mathieu Functions
-------------------------

 -- Function: int gsl_sf_mathieu_ce (int N, double Q, double X)
 -- Function: int gsl_sf_mathieu_ce_e (int N, double Q, double X, gsl_sf_result * RESULT)
 -- Function: int gsl_sf_mathieu_se (int N, double Q, double X)
 -- Function: int gsl_sf_mathieu_se_e (int N, double Q, double X, gsl_sf_result * RESULT)
     These routines compute the angular Mathieu functions ce_n(q,x) and se_n(q,x), respectively.

 -- Function: int gsl_sf_mathieu_ce_array (int NMIN, int NMAX, double Q, double X, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_se_array (int NMIN, int NMAX, double Q, double X, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
     These routines compute a series of the angular Mathieu functions ce_n(q,x) and se_n(q,x) of order n from NMIN to NMAX inclusive, storing the results in the array RESULT_ARRAY.

Radial Mathieu Functions
------------------------

 -- Function: int gsl_sf_mathieu_Mc (int J, int N, double Q, double X)
 -- Function: int gsl_sf_mathieu_Mc_e (int J, int N, double Q, double X, gsl_sf_result * RESULT)
 -- Function: int gsl_sf_mathieu_Ms (int J, int N, double Q, double X)
 -- Function: int gsl_sf_mathieu_Ms_e (int J, int N, double Q, double X, gsl_sf_result * RESULT)
     These routines compute the radial J-th kind Mathieu functions $Mc_n^(j)(q,x)$ Mc_n^{(j)}(q,x) and $Ms_n^(j)(q,x)$ Ms_n^{(j)}(q,x) of order N.

     The allowed values of J are 1 and 2. The functions for j = 3,4 can be computed as $M_n^(3) = M_n^(1) + iM_n^(2)$ M_n^{(3)} = M_n^{(1)} + iM_n^{(2)} and $M_n^(4) = M_n^(1) - iM_n^(2)$ M_n^{(4)} = M_n^{(1)} - iM_n^{(2)}, where $M_n^(j) = Mc_n^(j)$ M_n^{(j)} = Mc_n^{(j)} or $Ms_n^(j)$ Ms_n^{(j)}.

 -- Function: int gsl_sf_mathieu_Mc_array (int J, int NMIN, int NMAX, double Q, double X, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_Ms_array (int J, int NMIN, int NMAX, double Q, double X, gsl_sf_mathieu_workspace * WORK, double RESULT_ARRAY[])
     These routines compute a series of the radial Mathieu functions of kind J, with order from NMIN to NMAX inclusive, storing the results in the array RESULT_ARRAY.