The routines described in this section compute the Mathieu functions. The Mathieu functions are solutions of the two Mathieu's equations

     d^2x/d\phi^2 + (a - 2q\cos 2\phi)x = 0
     d^2x/du^2 - (a - 2q\cosh 2u)x = 0

Solutions to the first equation above are the angular Mathieu functions, whereas solutions to the second equation are the radial Mathieu functions.

The angular Mathieu functions can be divided into four types as

     x = \sum_(m=0)^\infty A_(2m+p) \cos(2m+p)\phi,   p = 0, 1,
     x = \sum_(m=0)^\infty B_(2m+p) \sin(2m+p)\phi,   p = 0, 1.

The nomenclature used for the angular Mathieu functions is ce_n for the first solution and se_n for the second.

Similar solutions exist for the radial Mathieu functions by replacing the trigonometric functions with their corresponding hyperbolic functions as shown below.

     x = \sum_(m=0)^\infty A_(2m+p) \cosh(2m+p)u,   p = 0, 1,
     x = \sum_(m=0)^\infty B_(2m+p) \sinh(2m+p)u,   p = 0, 1.

The nomenclature used for the radial Mathieu functions is Mc_n for the first solution and Ms_n for the second.  The hyperbolic series do not always converge at an acceptable rate.  Therefore most texts on the subject suggest using the following equivalent equations that are expanded in series of Bessel and Hankel functions.

     Mc_(2n)^(j)(x,q) = \sum_(m=0)^\infty (-1)^(r+k) A_(2m)^(2n)(q)
         [J_m(u_1)Z_m^(j)(u_2) + J_m(u_1)Z_m^(j)(u_2)]/A_2^(2n)
     Mc_(2n+1)^(j)(x,q) = \sum_(m=0)^\infty (-1)^(r+k) A_(2m+1)^(2n+1)(q)
         [J_m(u_1)Z_(m+1)^(j)(u_2) + J_(m+1)(u_1)Z_m^(j)(u_2)]/A_1^(2n+1)
     Ms_(2n)^(j)(x,q) = \sum_(m=1)^\infty (-1)^(r+k) B_(2m)^(2n)(q)
         [J_(m-1)(u_1)Z_(m+1)^(j)(u_2) + J_(m+1)(u_1)Z_(m-1)^(j)(u_2)]/B_2^(2n)
     Ms_(2n+1)^(j)(x,q) = \sum_(m=0)^\infty (-1)^(r+k) B_(2m+1)^(2n+1)(q)
         [J_m(u_1)Z_(m+1)^(j)(u_2) + J_(m+1)(u_1)Z_m^(j)(u_2)]/B_1^(2n+1)

where u_1 = \sqrtqe^-x and u_2 = \sqrtqe^x and

     Z_m^(1)(u) = J_m(u)
     Z_m^(2)(u) = Y_m(u)
     Z_m^(3)(u) = H_m^(1)(u)
     Z_m^(4)(u) = H_m^(2)(u)

where J_m(u), Y_m(u), H_m^(1)(u), and H_m^(2)(u) are the regular and irregular Bessel functions and the Hankel functions, respectively.

For more information, see Abromowitz and Stegun, Chapter 20.

The Mathieu functions are defined in the header file gsl_sf_mathieu.h.

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

 -- Function: int gsl_sf_mathieu_c (int N, double Q, double X, gsl_sf_result *RESULT)
 -- Function: int gsl_sf_mathieu_s (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: double *gsl_sf_mathieu_c_array (int NMIN, int NMAX, double Q, double X, gsl_sf_mathieu_workspace *WORK);
 -- Function: double *gsl_sf_mathieu_s_array (int NMIN, int NMAX, double QQ, double X, gsl_sf_mathieu_workspace *WORK);
     These routines compute a series of the angular Mathieu functions from NMIN to NMAX inclusive, storing the results in an array pointed to by the return pointer.  The calling routine is responsible for freeing the associated memory.

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

 -- Function: int gsl_sf_mathieu_mc_1 (int N, double Q, double X, gsl_sf_result *RESULT);
 -- Function: int gsl_sf_mathieu_ms_1 (int N, double Q, double X, gsl_sf_result *RESULT);
     These routines compute the radial Mathieu functions Mc_n^(1)(q,x) and Ms_n^(1)(q,x), respectively.

The remaining radial Mathieu functions (types 2, 3, and 4) have not yet been implemented.

Support Routines
----------------

 -- Function: gsl_sf_mathieu_workspace *gsl_sf_mathieu_alloc (const size_t N);
 -- Function: void gsl_sf_mathieu_free (gsl_sf_mathieu_workspace *WORKSPACE);
     These routines handle memory allocation and freeing for workspace arrays required by the ``_array'' versions of the Mathieu function routines.  The variable N specifies the maximum order of the functions to be computed, not the number of functions to be computed.