How to define Mc and Ms Mathieu functions in terms of MathieuC and MathieuS?

Reading multiple books about Mathieu functions, I always come across notation like $\operatorname{ce}_r(z,q)$, $\operatorname{se}_r(z,q)$ for angular functions and $\operatorname{Ce}_r(z,q)$, $\operatorname{Se}_r(z,q)$ as well as $\operatorname{Mc}_r^{(j)}(z,q)$, $\operatorname{Ms}_r^{(j)}(z,q)$ for radial functions. Example of the book with such notation is Abramovitz & Stegun "Handbook of mathematical functions".

For $\operatorname{ce}$ and $\operatorname{se}$ I figured that they can be defined in Mathematica as

ce[r_, z_, q_] = MathieuC[MathieuCharacteristicA[r, q], q, z];
se[r_, z_, q_] = MathieuS[MathieuCharacteristicB[r, q], q, z];

Similarly for $\operatorname{Ce}$ and $\operatorname{Se}$:

Ce[r_, z_, q_] = ce[r, I z, q];
Se[r_, z_, q_] = -I se[r, I z, q];

But how can I define $\operatorname{Mc}$ and $\operatorname{Ms}$?

Answer

Warning: long answer with little practical use.

As of Mathematica 10, there're no completely symbolic ways to define radial functions $\operatorname{Mc}^{(j)}$ and $\operatorname{Ms}^{(j)}$. Still, I'll show how to define them using as few numerics as I could achieve. In particular, here all the dependence of the functions on $z$ as in e.g. $\operatorname{Ms}^{(1)}(z,h)$ will be through built-in functions MathieuC and MathieuS.

I have to admit though, that although all this theoretically is valid, and should work, Mathematica appears to have a broken implementation of MathieuC and MathieuS. Thus, the radial functions for $z>2$ computed using these functions are garbage as of Mathematica 10.1 and earlier. It's still interesting for me to answer this question, so I proceed anyway.

First let's define the simplest Mathieu angular functions $\operatorname{ce}$ and $\operatorname{se}$. They'll be used in further definitions.

ce[r_, z_, q_] = MathieuC[MathieuCharacteristicA[r, q], q, z];
se[r_, z_, q_] = MathieuS[MathieuCharacteristicB[r, q], q, z];

These seem to be normalized already by the condition in DLMF 28.2.30. Mathematica can't evaluate these integrals symbolically, but I've checked them for multiple random parameters r and q to 60 digits to see that the condition is satisfied.

Next we'll want to define the "second solutions" of Mathieu equation, corresponding to the same characteristic value. The second solution corresponding to $\operatorname{ce}$ is $\operatorname{fe}$ and to $\operatorname{se}$ is $\operatorname{ge}$. Here we go for $\operatorname{fe}$ (see DLMF 28.5.1 and 28.5.5 for reference):

ϵ = 10^-14;
Qfe[n_, z_, q_] = MathieuS[MathieuCharacteristicA[n, q] (1 + ϵ), q, z];
QC[n_, q_] = (Qfe[n, 2 π, q] - Qfe[n, 0, q])/(2 π ce[n, 0, q]);
f[n_, z_, q_] = Qfe[n, z, q]/QC[n, q] - z ce[n, z, q];
ℭ[n_?NumericQ, q_?NumericQ] := Sqrt[π/NIntegrate[f[n, x, q]^2, {x, 0, 2 π}]]
fe[n_, z_, q_] := ℭ[n, q] (z ce[n, z, q] + f[n, z, q]);

Here Qfe denotes $$Q_n(q)\operatorname{fe}_n(z,q)=Q_n(q)\mathfrak{C}_n(q)(\operatorname{ce}_n(z,q)z+f_n(z,q)).$$

MMA's implementation of MathieuC and MathieuS functions appears such that at the values of MathieuCharacteristicA and MathieuCharacteristicB the opposite functions are "normalized" to 0. Thus we use ϵ to get close to those points, but not into them, and then restore the desired normalization. The $Q_n$ is the "unnormalization" factor, which we want to get rid of, and $\mathfrak{C}_n$ is what is denoted by $C_n$ in DLMF (because C is a used built-in MMA symbol). So we find this $C_n$ and use it to finally define the $\operatorname{fe}$. $C_n$ will be also used later on its own. To find $Q_n$ and $C_n$, we use $2\pi$-periodicity of $\operatorname{ce}_n$ and $f_n$ in $z$ and positivity of $\operatorname{ce}_n$ at $z=0$.

Analogously we do for $\operatorname{ge}$, using DLMF 28.5.2 and 28.5.5:

QgePrime[n_, z_, q_] = MathieuCPrime[MathieuCharacteristicB[n, q] (1 + ϵ), q, z];
Qge[n_, z_, q_] = MathieuC[MathieuCharacteristicB[n, q] (1 + ϵ), q, z];
QS[n_,q_]=(QgePrime[n,2π,q]-QgePrime[n,0,q])/(2π Derivative[0,1,0][se][n,0,q]);
g[n_, z_, q_] = Qge[n, z, q]/QS[n, q] - z se[n, z, q];
S[n_?NumericQ, q_?NumericQ] := Sqrt[π/NIntegrate[g[n, x, q]^2, {x, 0, 2 π}]]
ge[n_, z_, q_] := S[n, q] (z se[n, z, q] + g[n, z, q]);

Now, using DLMF 28.20(ii) we define simple solutions of modified Mathieu equation:

Ce[ν_, z_, q_] = ce[ν, I z, q];
Se[ν_, z_, q_] = -I se[ν, I z, q];

Fe[n_?IntegerQ, z_, q_] := -I fe[n, I z, q] /; n >= 0;
Ge[n_?IntegerQ, z_, q_] := ge[n, I z, q] /; n > 0;

Now we are close to start defining the radial Mathieu functions we were looking for in the first place. For this we have to find Fourier coefficients of $\operatorname{ce}$ and $\operatorname{se}$, i.e. $A_m^n(h)$ and $B_m^n(h)$ respectively. The expressions we'll use are consequences of DLMF 28.4(i); the explicit expressions can be seen in Abramovitz&Stegun 20.5.5.

A[m_?IntegerQ,n_?IntegerQ,h_?NumericQ]:=
  (1 - KroneckerDelta[0, m]/2) 1/π NIntegrate[ce[n, x, h] Cos[m x], {x, 0, 2 π}] /; 
  EvenQ[n] && EvenQ[m] || OddQ[n] && OddQ[m]
B[m_?IntegerQ,n_?IntegerQ,h_?NumericQ]:=1/π NIntegrate[se[n,x,h]Sin[m x],{x,0,2π}] /; 
  EvenQ[n] && EvenQ[m] || OddQ[n] && OddQ[m]

These coefficients won't be used to sum series though, it'd be too time consuming and prone to rounding errors. We'll use them as factors in normalization coefficients.

Now we need to define joining factors $g_{e,m}$ and $g_{o,m}$, as in DLMF 28.22.5-28.22.8:

gje[n_?EvenQ, h_] := (-1)^(n/2) Sqrt[2/π] ce[n, π/2, h^2]/A[0, n, h^2]
gje[n_?OddQ, h_] := (-1)^((n+1)/2) Sqrt[2/π] Derivative[0, 1, 0][ce][n, π/2, h^2]/(h A[1, n, h^2])
gjo[n_?OddQ, h_] := (-1)^((n-1)/2) Sqrt[2/π] se[n, π/2, h^2]/(h B[1, n, h^2])
gjo[n_?EvenQ, h_] := (-1)^(n/2) Sqrt[2/π] Derivative[0, 1, 0][se][n, π/2, h^2]/(h^2 B[2, n, h^2])

And, finally, we start defining the first types of radial functions, $\operatorname{Mc}^{(1)}$ and $\operatorname{Ms}^{(1)}$ via DLMF 28.22.1 and 28.22.2:

Mc1[m_,z_,h_]:=Sqrt[2/π] 1/(gje[m,h]ce[m, 0, h^2])Ce[m, z, h^2]
Ms1[m_,z_,h_]:=Sqrt[2/π] 1/(gjo[m,h]Derivative[0,1,0][se][m,0,h^2])Se[m,z,h^2]

The functions of second type, $\operatorname{Mc}^{(2)}$ and $\operatorname{Ms}^{(2)}$, are much harder to define. We'll use DLMF 28.22.3 and 28.22.4, but they together with joining factor equations 28.22.9 and 28.22.10 appear circular. One has to know the value of $\operatorname{Mc}^{(2)}_m(0,h)$ and, even harder to get, ${\operatorname{Ms}^{(2)}_m}'(0,h)$. To find the values, we'll use Abramovitz&Stegun 20.7.26. The related formulas in DLMF (28.28(ii)) appear to be broken as of this writing. So, for $\operatorname{Mc}^{(2)}_m(0,h)$:

Mc2z0[M_?OddQ, h_?NumericQ] := (-1)^((M-1)/2)/A[1,M,h^2](
  8 h)/π NIntegrate[(BesselY[1, h Sqrt[2 (1 + Cos[2 t])]] Cos[t])/(
    h Sqrt[2 (1 + Cos[2 t])]) ce[M, t, h^2], {t, 0, π/2}]
Mc2z0[M_?EvenQ, h_?NumericQ] := (-1)^(M/2)/A[0, M, h^2] 2/π NIntegrate[
   BesselY[0, h Sqrt[2 (1 + Cos[2 t])]] ce[M, t, h^2], {t, 0, π/2}]

Now the joining factor $f_{e,m}$ from DLMF 28.22.9:

fje[m_, h_] := -Sqrt[π/2] gje[m, h] Mc2z0[m, h]

And we are ready to define $\operatorname{Mc}^{(2)}$:

Mc2[m_, z_, h_] = Sqrt[2/π] 1/(gje[m, h] ce[m, 0, 
     h^2]) (2/(π ℭ[m, h^2]) Fe[m, z, h^2] - fje[m, h] Ce[m, z, h^2]);

Now $\operatorname{Ms}^{(2)}$ is harder, because DLMF 28.22.10 appears useless since integration in Abramovitz&Stegun 20.7.27 gives very imprecise numerical results in MMA near $z=0$, and switching integration and differentiation leads to apparently divergent (or at least hard to evaluate) integral. Thus we'll use DLMF 28.22.4, but with reference value of $\operatorname{Ms}^{(2)}$ computed via Abramovitz&Stegun 20.7.27 at $z=1/25$ (hopefully it doesn't vanish here for any given parameters $m$ and $h$). So, here're the reference integrals:

Ms2tA[M_?EvenQ, z_?NumericQ, h_?NumericQ] := (-1)^(M/2 + 1)/B[2, M, h^2] (
   8 h^2 Sinh[2 z])/π NIntegrate[(
   BesselY[2, h Sqrt[2 (Cosh[2 z] + Cos[2 t])]] Sin[2 t] se[M, t, 
     h^2])/(h^2 2 (Cosh[2 z] + Cos[2 t])), {t, 0, π/2}]

Ms2tA[M_?OddQ, z_?NumericQ, h_?NumericQ] := (-1)^((M - 1)/2)/B[1, M, h^2] (
   8 h Sinh[z])/π NIntegrate[(
   BesselY[1, h Sqrt[2 (Cosh[2 z] + Cos[2 t])]] Sin[t] se[M, t, 
     h^2])/(h Sqrt[2 (Cosh[2 z] + Cos[2 t])]), {t, 0, π/2}]

Then the joining factor $f_{o,m}$ from DLMF 28.22.10 is

fjo[m_, h_] := 
 With[{z = 1/25}, (Ms2tA[m, z, h] - Sqrt[2/π] 2/π (-Ge[m, z, h^2])/(
    gjo[m, h] Derivative[0, 1, 0][se][m, 0, h^2] S[m, h^2]))/-Ms1[m, z, h]]

Finally $\operatorname{Ms}^{(2)}$ can be defined:

Ms2[m_, z_, h_] = Sqrt[2/π] 1/(gjo[m, h] Derivative[0, 1, 0][se][m, 0, 
     h^2]) (-(2/(π S[m, h^2])) Ge[m, z, h^2] - fjo[m, h] Se[m, z, h^2]);