(* ::Package:: *)

BeginPackage["Mathieu`"]
Unprotect["Mathieu*"]
MathieuCE::usage="MathieuCE calculates the angular, first-kind, even solution to the Mathieu equation. 

MathieuCE[r,q,\[Theta]] takes integer values of r, real values of q, and real and imaginary values of \[Theta]."

MathieuCEPrime::usage="MathieuCEPrime calculates the first derivative of the angular, first-kind, even solution to the Mathieu equation, MathieuCE."

MathieuSE::usage="MathieuSE calculates the angular, first-kind, odd solution to the Mathieu equation.

MathieuSE[r,q,\[Theta]] takes positive integer values of r, real and imaginary q, and real or complex \[Theta]."

MathieuSEPrime::usage="MathieuSEPrime calculates the first derivative of the angular, first-kind, odd solution to the Mathieu equation, MathieuSE."

MathieuFE::usage="MathieuFE calculates the angular, second-kind, even solution to the Mathieu equation."

MathieuFEPrime::usage="MathieuFEPrime calculates the first derivative of the angular, second-kind, even solution to the Mathieu equation, MathieuFE."

MathieuFO::usage="MathieuFO calculates the angular, second-kind, odd solution to the Mathieu equation."

MathieuFOPrime::usage="MathieuFOPrime calculates the first derivative of the angular, second-kind, odd solution to the Mathieu equation, MathieuFO."

MathieuJE::usage="MathieuJE calculates the radial, first-kind, even solution to the modified Mathieu equation.

MathieuJE takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuJEPrime::usage="MathieuJEPrime calculates the first TEHE derivatve of the radial, first-kind, even solution to the modified Mathieu equation, MathieuJE."

MathieuJO::usage="MathieuJO calculates the radial, first-kind, odd solution to the modified Mathieu equation.

MathieuJO takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuJOPrime::usage="MathieuJOPrime calculates the first derivatve of the radial, first-kind, odd solution to the modified Mathieu equation, MathieuJO."

MathieuIE::usage="MathieuIE calculates the radial, first-kind, even, modified solution to the modified Mathieu equation, analogous to BesselI.

MathieuJE takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuIEPrime::usage="MathieuIEPrime calculates the first derivatve of the radial, first-kind, even solution to the modified Mathieu equation, MathieuIE."

MathieuIO::usage="MathieuIO calculates the radial, first-kind, odd, modified solution to the modified Mathieu equation, analogous to BesselI.

MathieuIO takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuIOPrime::usage="MathieuIOPrime calculates the first derivatve of the radial, first-kind, odd solution to the modified Mathieu equation, MathieuIO."

MathieuYE::usage="MathieuYE calculates the radial, second-kind, even solution to the modified Mathieu equation.

MathieuYE takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuYEPrime::usage="MathieuYEPrime calculates the first derivatve of the radial, second-kind, even solution to the modified Mathieu equation, MathieuYE."

MathieuYO::usage="MathieuYO calculates the radial, second-kind, odd solution to the modified Mathieu equation.

MathieuYO takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuYOPrime::usage="MathieuYOPrime calculates the first derivatve of the radial, second-kind, odd solution to the modified Mathieu equation, MathieuYO."

MathieuKE::usage="MathieuKE calculates the radial, third-kind, even solution to the modified Mathieu equation.

MathieuKE takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuKEPrime::usage="MathieuKEPrime calculates the first derivatve of the radial, third-kind, even solution to the modified Mathieu equation, MathieuKE."

MathieuKO::usage="MathieuKO calculates the radial, third-kind, odd solution to the modified Mathieu equation.

MathieuKO takes positive, integer values of r, real q, and real or complex \[Theta]."

MathieuKOPrime::usage="MathieuKOPrime calculates the first derivatve of the radial, third-kind, odd solution to the modified Mathieu equation, MathieuKO."

MathieuAngularEvenWronskian::usage="MathieuAngularEvenWronskian[r,q,\[Theta]] calculates the Wronskian of the even angular solutions to the Mathieu equation, which should be 2/\[Pi]."

MathieuAngularOddWronskian::usage="MathieuAngularOddWronskian[r,q,\[Theta]] calculates the Wronskian of the odd angular solutions to the Mathieu equation, which should be 2/\[Pi]."

MathieuRadialEvenWronskian::usage="MathieuRadialEvenWronskian[r,q,\[Theta]] calculates the Wronskian of the even radial solutions to the Mathieu equation, which should be 2/\[Pi]."

MathieuRadialOddWronskian::usage="MathieuRadialOddWronskian[r,q,\[Theta]] calculates the Wronskian of the odd radial solutions to the Mathieu equation, which should be 2/\[Pi]."

MathieuModifiedEvenWronskian::usage="MathieuModifiedEvenWronskian[r,q,\[Theta]] calculates the Wronskian of the even radial solutions to the Mathieu equation, which should be 1."

MathieuModifiedOddWronskian::usage="MathieuModifiedOddWronskian[r,q,\[Theta]] calculates the Wronskian of the odd radial solutions to the Mathieu equation, which should be 1."


MathieuClearCache::usage="Clears internal caches to save memory."


Begin["`Private`"]
 
isZero[x_]:=x==0
OddZero[x_]:=OddQ[x]\[Or]x==0
EvenNotZero[x_]:=EvenQ[x]\[And]x!=0
$RecursionLimit=10000

isInteger[x_]:=x==Round[x]

(*creates mathieu and error generator function*)

makeMathieuFunction[name_]:=
(MakeExpression[RowBox[{StringJoin["Mathieu",ToString[name]],"[r_?isInteger,q_,\[Theta]_]",":=",ToString[name],"[Round[r],SetPrecision[q,Min[Precision[q],Precision[\[Theta]]]],\[Theta]]","/;","Min[Precision[q],Precision[\[Theta]]]","<","Infinity"}],StandardForm])

MathieuSE[r_?isZero,q_,\[Theta]_]:=Message[MathieuSE::zord]
MathieuFO[r_?isZero,q_,\[Theta]_]:=Message[MathieuFO::zord]
MathieuJO[r_?isZero,q_,\[Theta]_]:=Message[MathieuJO::zord]
MathieuYO[r_?isZero,q_,\[Theta]_]:=Message[MathieuYO::zord]
MathieuIO[r_?isZero,q_,\[Theta]_]:=Message[MathieuIO::zord]
MathieuKO[r_?isZero,q_,\[Theta]_]:=Message[MathieuKO::zord]
MathieuSEPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuSE::zord]
MathieuFOPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuFO::zord]
MathieuJOPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuJO::zord]
MathieuYOPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuYO::zord]
MathieuIOPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuIO::zord]
MathieuKOPrime[r_?isZero,q_,\[Theta]_]:=Message[MathieuKO::zord]

MathieuSE::zord= "There is no zero-order MathieuSE . "
MathieuFO::zord= "There is no zero-order MathieuFO . "
MathieuJO::zord= "There is no zero-order MathieuJO . "
MathieuYO::zord= "There is no zero-order MathieuYO . "
MathieuIO::zord= "There is no zero-order MathieuIO . "
MathieuKO::zord= "There is no zero-order MathieuKO . "
shorthandFN:={CE, SE, FE, FO, JE, JO, YE, YO, IE, IO, KE, KO, CEPrime, SEPrime, FEPrime, FOPrime,JEPrime,
JOPrime, YEPrime, YOPrime, IEPrime, IOPrime, KEPrime, KOPrime}


(*creates mathieu top level function*)
ReleaseHold[(makeMathieuFunction/@shorthandFN)];
MathieuAngularEvenWronskian[r_,q_,\[Theta]_]:=MathieuCE[r,q,\[Theta]] MathieuFEPrime[r,q,\[Theta]]-MathieuCEPrime[r,q,\[Theta]] MathieuFE[r,q,\[Theta]]
MathieuAngularOddWronskian[r_,q_,\[Theta]_]:=MathieuSEPrime[r,q,\[Theta]] MathieuFO[r,q,\[Theta]]-MathieuSE[r,q,\[Theta]] MathieuFOPrime[r,q,\[Theta]]
MathieuRadialEvenWronskian[r_,q_,\[Theta]_]:=MathieuJE[r,q,\[Theta]] MathieuYEPrime[r,q,\[Theta]]-MathieuJEPrime[r,q,\[Theta]] MathieuYE[r,q,\[Theta]]
MathieuRadialOddWronskian[r_,q_,\[Theta]_]:=MathieuJO[r,q,\[Theta]] MathieuYOPrime[r,q,\[Theta]]-MathieuJOPrime[r,q,\[Theta]] MathieuYO[r,q,\[Theta]]
MathieuModifiedEvenWronskian[r_,q_,\[Theta]_]:=MathieuIEPrime[r,q,\[Theta]] MathieuKE[r,q,\[Theta]]-MathieuIE[r,q,\[Theta]] MathieuKEPrime[r,q,\[Theta]]
MathieuModifiedOddWronskian[r_,q_,\[Theta]_]:=MathieuIOPrime[r,q,\[Theta]] MathieuKO[r,q,\[Theta]]-MathieuIO[r,q,\[Theta]] MathieuKOPrime[r,q,\[Theta]]


(*efficiency functions Angular Part*)
AnPNCoE[r_] := If[r<0,-1,1]
AnPNCoER[r_,q_] := If[Re[q]<0,(-1)^(r/2),1]
AnPNCoERR[r_,q_] := If[Re[q]<0,(-1)^((r-1)/2),1]
AnPNCoERRR[r_,q_] := If[Re[q]<0,(-1)^(r/2+1),1]

(*efficiency functions Radial Part*)
RaNeg[r_] := If[Negative[r],(-1)^r,1]
RaPNCoE[r_] := (-1)^(Abs[r]/2)
RaPNCoER[r_] := (-1)^((Abs[r]-1)/2)
RaHOCoE[r_] := If[r==0,1/2,1]


(*With flexible limits*)
$MathieuLowerLimit:=20
lowerLimitFunction[l_,q_] := 0
upperLimitFunction[l_,q_] := Max[$MathieuLowerLimit, Round[l+2*Sqrt[Abs[q]]]]
mySum[element_,arg_,ArgL_,ArgQ_]:= (Sum[ReleaseHold[element],{arg,lowerLimitFunction[ArgL,ArgQ],upperLimitFunction[ArgL,ArgQ]}])

(*Coefficients*)

(*General coefficients*)
characteristicValue[r_,q_,evenodd_,j_,c_]:=If[c==0,If[evenodd==0,MathieuCharacteristicA[r,q],MathieuCharacteristicB[r,q]],If[evenodd==0,a[r,q,j+1],b[r,q,j+1]]]
genCoefficient[n_,r_,q_,evenodd_,p_,j_,c_]:=genCoefficient[n,r,q,evenodd,p,j,c]=If[n<=r,(characteristicValue[r,q,evenodd,j,c]-(n)^2)/q-1/genCoefficient[n-2,r,q,evenodd,p,j,c],-q/((n+2)^2-characteristicValue[r,q,evenodd,j,c]+q genCoefficient[n+2,r,q,evenodd,p,j,c])]

(*First-Kind coefficients*)
(*The last two parameters are for the series expansion; for numerical evaluation, c parameter is zero and j is unused*)
A[m_,r_?EvenQ,q_,j_,c_]:=Ae[2m,r,q,j,c]
Ae[0,r_,q_,j_,c_]:=1
Ae[n_,r_,q_,j_,c_]:=Ae[n,r,q,j,c]=Ae[n-2,r,q,j,c]*genCoefficient[n-2,r,q,0,0,j,c]

A[m_,r_?OddQ,q_,j_,c_]:=Ao[2m+1,r,q,j,c]
Ao[1,r_,q_,j_,c_]:=1
Ao[n_,r_,q_,j_,c_]:=Ao[n,r,q,j,c]=Ao[n-2,r,q,j,c]*genCoefficient[n-2,r,q,0,1,j,c]

(*Ratio between first-kind,even coefficients*)
genCoefficient[0,r_,q_,0,0,j_,c_]:=characteristicValue[r,q,0,j,c]/q
genCoefficient[2,r_,q_,0,0,j_,c_]:=(characteristicValue[r,q,0,j,c]-4)/q-2/genCoefficient[0,r,q,0,0,j,c]
genCoefficient[x_,r_,q_,0,0,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]
genCoefficient[1,r_,q_,0,1,j_,c_]:=(characteristicValue[r,q,0,j,c]-1)/q-1
genCoefficient[3,r_,q_,0,1,j_,c_]:=(characteristicValue[r,q,0,j,c]-9)/q-1/genCoefficient[1,r,q,0,1,j,c]
genCoefficient[x_,r_,q_,0,1,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]+1

(*First-Kind,Odd Coefficients*)
B[m_,r_?EvenQ,q_,j_,c_]:=Be[2*m,r,q,j,c]
Be[0,r_,q_,j_,c_]:=0
Be[2,r_,q_,j_,c_]:=1
Be[n_,r_,q_,j_,c_]:=Be[n,r,q,j,c]=Be[n-2,r,q,j,c]*genCoefficient[n-2,r,q,1,0,j,c]

B[m_,r_?OddQ,q_,j_,c_]:=Bo[2*m+1,r,q,j,c]
Bo[1,r_,q_,j_,c_]:=1
Bo[n_,r_,q_,j_,c_]:=Bo[n,r,q,j,c]=Bo[n-2,r,q,j,c]*genCoefficient[n-2,r,q,1,1,j,c]

(*Ratio between first-kind,odd coefficients*)
genCoefficient[0,r_,q_,1,0,j_,c_]:=0
genCoefficient[2,r_,q_,1,0,j_,c_]:=(characteristicValue[r,q,1,j,c]-4)/q
genCoefficient[x_,r_,q_,1,0,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]
genCoefficient[1,r_,q_,1,1,j_,c_]:=1+(characteristicValue[r,q,1,j,c]-1)/q
genCoefficient[3,r_,q_,1,1,j_,c_]:=(characteristicValue[r,q,1,j,c]-9)/q-1/genCoefficient[1,r,q,1,1,j,c]
genCoefficient[x_,r_,q_,1,1,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]+1

(*Second-Kind,Even coefficients*)
Gnorm[m_,r_,q_,j_,c_]:=Gnorm[m,r,q,j,c]=G[m,r,q,j,c]/(mySum[Hold[A[\[Nu],r,q,j,c]],\[Nu],Abs[r]/2,q])

(*Un-normalized coefficients*)
G[m_,r_?EvenQ,q_,j_,c_]:=G[m,r?EvenQ,q,j,c]=Ge[2m,r,q,j,c]
G[m_,r_?OddQ,q_,j_,c_]:=G[m,r?OddQ,q,j,c]=Go[2m+1,r,q,j,c]

Ge[0,r_,q_,j_,c_]:=0
Ge[n_,r_,q_,j_,c_]:=(Qe[n,r,q,j,c]-\[Rho]qe[n,r,q,j,c] Ae[n,r,q,j,c])
Go[n_,r_,q_,j_,c_]:=(Qo[n,r,q,j,c]-\[Rho]qo[n,r,q,j,c] Ao[n,r,q,j,c])

(*Scalar to normalize the first-kind coefficients*)
\[Rho]qe[n_,r_,q_,j_,c_]:=(((1/(2Ae[0,r,q,j,c]))((If[c==0,MathieuCharacteristicA[r,q],a[r,q,j+1]]-4) Qe[2,r,q,j,c]/q-Qe[4,r,q,j,c]))-2If[c==0,MathieuCharacteristicA[r,q],a[r,q,j+1]]/(q^2))
\[Rho]qo[n_,r_,q_,j_,c_]:=((1/(2Ao[1,r,q,j,c]))((If[c==0,MathieuCharacteristicA[r,q],a[r,q,j+1]]-1+q) Qo[1,r,q,j,c]/q-Qo[3,r,q,j,c])-(1/q))

(*Qe,Qo are ratios between G,analogous to Ve,Vo and We,Wo*)
Qe[x_,r_,q_,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]
Qe[x_,r_,q_,j_,c_]:=-2 (x+2) Ae[(x+2),r,q,j,c]/(q)/;x>=2*upperLimitFunction[Abs[r]/2,q]-2
Qe[n_,r_,q_,j_,c_]:=Qe[n,r,q,j,c]=(((If[c==0,MathieuCharacteristicA[r,q],a[r,q,j+1]]-(n+2)^2) Qe[n+2,r,q,j,c]-2 (n+2) Ae[n+2,r,q,j,c])/q-Qe[n+4,r,q,j,c])
Qo[x_,r_,q_,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]+1
Qo[x_,r_,q_,j_,c_]:=-2(x+2) Ao[x+2,r,q,j,c]/q/;x>=2*upperLimitFunction[Abs[r]/2,q]-1
Qo[n_,r_,q_,j_,c_]:=Qo[n,r,q,j,c]=(((If[c==0,MathieuCharacteristicA[r,q],a[r,q,j+1]]-(n+2)^2) Qo[n+2,r,q,j,c]-2 (n+2) Ao[n+2,r,q,j,c])/q-Qo[n+4,r,q,j,c])

(*Second-Kind,odd coefficients*)
Hnorm[m_,r_?EvenQ,q_,j_,c_]:=H[m,r,q,j,c]/(mySum[Hold[2 \[Nu] B[\[Nu],r,q,j,c]],\[Nu],Abs[r]/2,q])
Hnorm[m_,r_?OddQ,q_,j_,c_]:=H[m,r,q,j,c]/(mySum[Hold[(2 \[Nu]+1)B[\[Nu],r,q,j,c]],\[Nu],Abs[r]/2,q])

(*Un-normalized coefficients*)
H[m_,r_?EvenQ,q_,j_,c_]:=He[2m,r,q,j,c]
H[m_,r_?OddQ,q_,j_,c_]:=Ho[2m+1,r,q,j,c]

He[n_,r_,q_,j_,c_]:=(Te[n,r,q,j,c]-\[Rho]te[n,r,q,j,c] Be[n,r,q,j,c])
Ho[n_,r_,q_,j_,c_]:=(To[n,r,q,j,c]-\[Rho]to[n,r,q,j,c] Bo[n,r,q,j,c])

(*Scalar to normalize the first-kind coefficients*)
\[Rho]te[n_,r_,q_,j_,c_]:=(1/Be[2,r,q,j,c])(Te[2,r,q,j,c]-If[c==0,MathieuCharacteristicB[r,q],b[r,q,j+1]]Te[0,r,q,j,c]/q)
\[Rho]to[n_,r_,q_,j_,c_]:=(1/(2 Bo[1,r,q,j,c]))(To[3,r,q,j,c]-(If[c==0,MathieuCharacteristicB[r,q],b[r,q,j+1]]-1-q)To[1,r,q,j,c]/q)-(1/q)

Te[x_,r_,q_,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]
Te[0,r_,q_,j_,c_]:=(4Be[2,r,q,j,c]+(If[c==0,MathieuCharacteristicB[r,q],b[r,q,j+1]]-4)Te[2,r,q,j,c])/(2 q)-Te[4,r,q,j,c]/2
Te[x_,r_,q_,j_,c_]:=2 (x+2) Be[(x+2),r,q,j,c]/q/;x>=2*upperLimitFunction[Abs[r]/2,q]-2
Te[n_,r_,q_,j_,c_]:=Te[n,r,q,j,c]=((If[c==0,MathieuCharacteristicB[r,q],b[r,q,j+1]]-(n+2)^2) Te[n+2,r,q,j,c]+2 (n+2) Be[n+2,r,q,j,c])/q-Te[n+4,r,q,j,c](*/2*)
To[x_,r_,q_,j_,c_]:=0/;x>=2*upperLimitFunction[Abs[r]/2,q]
To[x_,r_,q_,j_,c_]:=2(x+2) Bo[x+2,r,q,j,c]/q/;x>=2upperLimitFunction[Abs[r]/2,q]-1
To[n_,r_,q_,j_,c_]:=To[n,r,q,j,c]=((If[c==0,MathieuCharacteristicB[r,q],b[r,q,j+1]]-(n+2)^2) To[n+2,r,q,j,c]+2 (n+2) Bo[n+2,r,q,j,c])/q-To[n+4,r,q,j,c]


(*Angular Functions*)

(*Normalization Constant*)
NormalizationSum[r_,q_,AB_,first_,squared_,squareroot_,prefactor_]:=NormalizationSum[r,q,AB,first,squared,squareroot,prefactor]=(mySum[Hold[(prefactor)*(AB[m,r,q,j,0])^squared],m,Abs[r]/2,q]+If[first==1,AB[0,r,q,j,0],0])^squareroot


cutoff = 10^-5

(*Angular,first-kind,even Mathieu function*)
JFEven[r_,q_]:=JFEven[r,q]=CE[r,q,0]/JE[r,q,0]
CE[r_,q_?isZero,\[Theta]_]:=Cos[r \[Theta]]
CE[r_?EvenQ,q_,\[Theta]_?NumberQ]:=AnPNCoER[r,q] mySum[Hold[A[m,r,q,j,0]Cos[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,A,1,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
CE[r_?OddQ,q_,\[Theta]_?NumberQ]:=AnPNCoERR[r,q]mySum[Hold[A[m,r,q,j,0]Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,A,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
CE[1,q_,\[Theta]_?NumberQ]:=JE[1,q,-I \[Theta]]/JE[1,q,0]/;(Im[\[Theta]]==0)&&(Abs[Re[q]]<cutoff)
CE[r_,q_,\[Theta]_?NumberQ]:=JE[r,q,-I \[Theta]](JFEven[r,q])/;Im[\[Theta]]!=0

(*Angular,first-kind,odd Mathieu function*)
JFOdd[r_,q_]:=JFOdd[r,q]=SEPrime[r,q,0]/JOPrime[r,q,0]
SE[r_,q_?isZero,\[Theta]_]:=Sin[r \[Theta]]
SE[r_?EvenQ,q_,\[Theta]_]:=AnPNCoE[r]AnPNCoERRR[r,q]mySum[Hold[B[m,r,q,j,0]Sin[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,B,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
SE[r_?OddQ,q_,\[Theta]_]:=AnPNCoE[r]AnPNCoERR[r,q]mySum[Hold[B[m,r,q,j,0]Sin[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,B,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
SE[1,q_,\[Theta]_]:=AnPNCoE[r]I JO[1,q,-I \[Theta]]/(JOPrime[1,q,0])/;(Im[\[Theta]]==0)&&(Abs[Re[q]]<cutoff)
SE[r_,q_,\[Theta]_]:=AnPNCoE[r]I JO[r,q,-I \[Theta]](JFOdd[r,q])/;Im[\[Theta]]!=0

(*First derivative of the angular,first-kind,even Mathieu function*)
CEPrime[r_,z_?isZero,\[Theta]_]:=-r Sin[r \[Theta]]
CEPrime[r_?EvenQ,q_,\[Theta]_?NumberQ]:=AnPNCoER[r,q]mySum[Hold[A[m,r,q,j,0](2m)(-Sin[2m \[Theta]])],m,Abs[r]/2,q]/NormalizationSum[r,q,A,1,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
CEPrime[r_?OddQ,q_,\[Theta]_?NumberQ]:=AnPNCoERR[r,q]mySum[Hold[A[m,r,q,j,0](2m+1)(-Sin[(2m+1)\[Theta]])],m,Abs[r+1]/2,q]/NormalizationSum[r,q,A,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
CEPrime[1,q_,\[Theta]_?NumberQ]:=-I JEPrime[1,q,-I \[Theta]]/JE[1,q,0]/;(Im[\[Theta]]==0)&&(Abs[Re[q]]<cutoff)
CEPrime[r_,q_,\[Theta]_?NumberQ]:=-I JEPrime[r,q,\[Theta]/I]JFEven[r,q]/;Im[\[Theta]]!=0

(*First derivative of the angular,first-kind,odd Mathieu function*)
SEPrime[r_,q_?isZero,\[Theta]_]:=r Cos[r \[Theta]]
SEPrime[r_?EvenQ,q_,\[Theta]_]:=AnPNCoERRR[r,q]mySum[Hold[B[m,r,q,j,0](2m)Cos[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,B,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
SEPrime[r_?OddQ,q_,\[Theta]_]:=AnPNCoERR[r,q]mySum[Hold[B[m,r,q,j,0](2m+1)Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,B,0,2,1/2,1]/;Im[\[Theta]]==0&&(r!=1 || Abs[Re[q]]>=cutoff)
SEPrime[1,q_,\[Theta]_]:=JOPrime[1,q,-I \[Theta]]/(JOPrime[1,q,0])/;(Im[\[Theta]]==0)&&(Abs[Re[q]]<cutoff)
SEPrime[r_,q_,\[Theta]_]:=JOPrime[r,q,-I \[Theta]] JFOdd[r,q]/;Im[\[Theta]]!=0

(*Angular,second-kind,even Mathieu function*)
FE[r_,q_?isZero,\[Theta]_]:=2/\[Pi] Sin[r \[Theta]]/r
FE[r_?EvenQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoER[r,q](\[Theta] CEnorm[r,q,\[Theta]]+mySum[Hold[Gnorm[m,r,q,j,0]Sin[2m \[Theta]]],m,Abs[r]/2,q])NormalizationSum[r,q,A,1,2,1/2,1]/(NormalizationSum[r,q,A,0,1,1,1](1+NormalizationSum[r,q,Gnorm,0,1,1,2m]))
FE[r_?OddQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERR[r,q](\[Theta] CEnorm[r,q,\[Theta]]+mySum[Hold[Gnorm[m,r,q,j,0]Sin[(2m+1)\[Theta]]],m,Abs[r+1]/2,q])NormalizationSum[r,q,A,0,2,1/2,1]/((1+NormalizationSum[r,q,Gnorm,0,1,1,2m+1])(NormalizationSum[r,q,A,0,1,1,1]))
(*FE[r_,q_,\[Theta]_?NumberQ]:=YE[r,q,-I \[Theta]](FE[r,q,0]/YE[r,q,0])/;Im[\[Theta]]\[NotEqual]0*)

(*Angular,second-kind,odd Mathieu function*)
FO[r_,q_?isZero,\[Theta]_]:=2/\[Pi] Cos[r \[Theta]]/r
FO[r_?EvenQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERRR[r,q](\[Theta] SEnorm[r,q,\[Theta]]+mySum[Hold[Hnorm[m,r,q,j,0]Cos[2m \[Theta]]],m,Abs[r]/2,q])NormalizationSum[r,q,B,0,2,1/2,1]/((NormalizationSum[r,q,B,0,1,1,2m])(NormalizationSum[r,q,Hnorm,0,1,1,1]))
FO[r_?OddQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERR[r,q](\[Theta] SEnorm[r,q,\[Theta]]+mySum[Hold[Hnorm[m,r,q,j,0]Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q])NormalizationSum[r,q,B,0,2,1/2,1]/((NormalizationSum[r,q,Hnorm,0,1,1,1])(NormalizationSum[r,q,B,0,1,1,2m+1]))

(*First derivative of the angular,second-kind,even Mathieu function*)
FEPrime[r_?EvenQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoER[r,q](CEnorm[r,q,\[Theta]]+\[Theta] CEPrimenorm[r,q,\[Theta]]+mySum[Hold[2m Gnorm[m,r,q,j,0]Cos[2m \[Theta]]],m,Abs[r]/2,q])NormalizationSum[r,q,A,1,2,1/2,1]/(NormalizationSum[r,q,A,0,1,1,1](1+NormalizationSum[r,q,Gnorm,0,1,1,2m]))
FEPrime[r_?OddQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERR[r,q](CEnorm[r,q,\[Theta]]+\[Theta] CEPrimenorm[r,q,\[Theta]]+mySum[Hold[(2m+1)Gnorm[m,r,q,j,0]Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q])NormalizationSum[r,q,A,0,2,1/2,1]/((1+NormalizationSum[r,q,Gnorm,0,1,1,2m+1])(NormalizationSum[r,q,A,0,1,1,1]))

(*First derivative of the angular,second-kind,odd Mathieu function*)
FOPrime[r_?EvenQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERRR[r,q](\[Theta] SEPrimenorm[r,q,\[Theta]]+SEnorm[r,q,\[Theta]]+mySum[Hold[Hnorm[m,r,q,j,0](-2m)Sin[2m \[Theta]]],m,Abs[r]/2,q])NormalizationSum[r,q,B,0,2,1/2,1]/((NormalizationSum[r,q,B,0,1,1,2m])(NormalizationSum[r,q,Hnorm,0,1,1,1]))
FOPrime[r_?OddQ,q_,\[Theta]_?NumberQ]:=2/\[Pi] AnPNCoERR[r,q](\[Theta] SEPrimenorm[r,q,\[Theta]]+SEnorm[r,q,\[Theta]]+mySum[Hold[Hnorm[m,r,q,j,0](2m+1)(-Sin[(2m+1)\[Theta]])],m,Abs[r+1]/2,q])NormalizationSum[r,q,B,0,2,1/2,1]/((NormalizationSum[r,q,Hnorm,0,1,1,1])(NormalizationSum[r,q,B,0,1,1,2m+1]))

(*re-normalized CE,SE,CEPrime,and SEPrime*)
CEnorm[r_?EvenQ,q_,\[Theta]_?NumberQ]:=mySum[Hold[A[m,r,q,j,0]Cos[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,A,0,1,1,1]
CEnorm[r_?OddQ,q_,\[Theta]_?NumberQ]:=mySum[Hold[A[m,r,q,j,0]Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,A,0,1,1,1]
SEnorm[r_?EvenQ,q_,\[Theta]_]:=mySum[Hold[B[m,r,q,j,0]Sin[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,B,0,1,1,2m]
SEnorm[r_?OddQ,q_,\[Theta]_]:=mySum[Hold[B[m,r,q,j,0]Sin[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,B,0,1,1,2m+1]

CEPrimenorm[r_?EvenQ,q_,\[Theta]_?NumberQ]:=mySum[Hold[A[m,r,q,j,0](2m)(-Sin[2m \[Theta]])],m,Abs[r]/2,q]/NormalizationSum[r,q,A,0,1,1,1]
CEPrimenorm[r_?OddQ,q_,\[Theta]_?NumberQ]:=mySum[Hold[A[m,r,q,j,0](2m+1)(-Sin[(2m+1)\[Theta]])],m,Abs[r+1]/2,q]/NormalizationSum[r,q,A,0,1,1,1]
SEPrimenorm[r_?EvenQ,q_,\[Theta]_]:=mySum[Hold[B[m,r,q,j,0](2m)Cos[2m \[Theta]]],m,Abs[r]/2,q]/NormalizationSum[r,q,B,0,1,1,2m]
SEPrimenorm[r_?OddQ,q_,\[Theta]_]:=mySum[Hold[B[m,r,q,j,0](2m+1)Cos[(2m+1)\[Theta]]],m,Abs[r+1]/2,q]/NormalizationSum[r,q,B,0,1,1,2m+1]


(*Radial Functions*)

(*General recursion to speed up calculation of the standard Bessel functions*)
makeBesselTable[besselfunc_,tablefunc_,r_?isInteger,z_,qsign_,Q_]:=
(tablefunc[0]=besselfunc[0,z];
tablefunc[1]=besselfunc[1,z];
tablefunc[n_]:=If[n<0,(-qsign)^n tablefunc[-n],tablefunc[n]=-qsign tablefunc[n-2]+(2 (n-1)/z) tablefunc[n-1]])/;((besselfunc===BesselY||besselfunc===BesselK)&&Im[z]==0&&z!=0)

makeBesselTable[besselfunc_,tablefunc_,r_?isInteger,z_,qsign_,Q_]:=
(z2=SetPrecision[z,Precision[z]+30];
tablefunc[upperLimitFunction[Abs[r]/2,Q]+r]=besselfunc[upperLimitFunction[Abs[r]/2,Q]+r,z2];
tablefunc[upperLimitFunction[Abs[r]/2,Q]+r-1]=besselfunc[upperLimitFunction[Abs[r]/2,Q]+r-1,z2];
x=2; While[besselfunc[upperLimitFunction[Abs[r]/2,Q]+r-x,z2]==0,tablefunc[upperLimitFunction[Abs[r]/2,Q]+r-x-1]=besselfunc[upperLimitFunction[Abs[r]/2,Q]+r-x-1,z2];x+=1];
tablefunc[n_]:=If[n<0,(-qsign)^n tablefunc[-n],tablefunc[n]=-qsign tablefunc[n+2]+(2 (n+1)/z) tablefunc[n+1]])/;(besselfunc===BesselJ||besselfunc===BesselI&&Im[z]==0&&z!=0)

makeBesselTable[besselfunc_,tablefunc_,r_?isInteger,z_,qsign_,Q_]:=
(tablefunc[n_]:=besselfunc[n,z])/;Im[z]!=0||z==0

(*General form of the radial Mathieu functions,using temporary variables to calculate the needed Bessel functions*)

(*Work around bugs in Mathematica versions after 7.0*)
imagTruncate[z_]:=If[Im[z]==0,Re[z],z]

RadialE[r_,q_,z_,bessela_,besselb_,modifieda_,modifiedb_,l_,p_,qsign_]:=
Block[
{sqrtq=Sqrt[qsign*q],emu=Exp[imagTruncate[z]],functiona,functionb,A},
makeBesselTable[bessela,functiona,r,sqrtq/emu,qsign,q];
makeBesselTable[besselb,functionb,r,sqrtq*emu,qsign,q];
(1/(A[l,r,q,j,0]))mySum[Hold[(modifieda A[m,r,q,j,0])(functiona[m-l]functionb[m+l+p]+modifiedb functiona[m+l+p] functionb[m-l])],m,l,q]
]

RadialO[r_,q_,z_,bessela_,besselb_,modifieda_,modifiedb_,l_,p_,qsign_]:=
Block[
{sqrtq=Sqrt[qsign*q],emu=Exp[imagTruncate[z]],functiona,functionb,B},
makeBesselTable[bessela,functiona,r,sqrtq/emu,qsign,q];
makeBesselTable[besselb,functionb,r,sqrtq*emu,qsign,q];
(1/(B[l,r,q,j,0])) mySum[Hold[modifieda B[m,r,q,j,0](functiona[m-l] functionb[m+l+p]-modifiedb*functiona[m+l+p] functionb[m-l])],m,l,q]
]

(*General form of the derivatives of Bessel functions,using our recursion again for speed*)
besselprime[BesselK,function_,qsign_]:=
Function[n,-(function[n-1]-qsign function[n+1])/2]
besselprime[bessel_,function_,qsign_]:=
Function[n,(function[n-1]-qsign function[n+1])/2]

(*General form of the derivatives of the radial Mathieu functions*)
RadialEPrime[r_,q_,z_,bessela_,besselb__,modifieda_,modifiedb_,l_,p_,qsign_]:=
Block[
{sqrtq=Sqrt[qsign*q],emu=Exp[imagTruncate[z]],functiona,functionb,functionaprime,functionbprime,A},
makeBesselTable[bessela,functiona,r+1,sqrtq/emu,qsign,q];
makeBesselTable[besselb,functionb,r+1,sqrtq*emu,qsign,q];
functionaprime=besselprime[bessela,functiona,qsign];
functionbprime=besselprime[besselb,functionb,qsign];
(1/(A[l,r,q,j,0])) mySum[Hold[modifieda*A[m,r,q,j,0] ((sqrtq*emu*functiona[m-l] functionbprime[m+l+p])-(sqrtq/emu functionb[m+l+p] functionaprime[m-l])+modifiedb*((sqrtq emu functiona[m+l+p] functionbprime[m-l])-(sqrtq/emu functionb[m-l] functionaprime[m+l+p])))],m,l,q]
]

RadialOPrime[r_,q_,z_,bessela_,besselb__,modifieda_,modifiedb_,l_,p_,qsign_]:=
Block[
{sqrtq=Sqrt[qsign*q],emu=Exp[imagTruncate[z]],functiona,functionb,functionaprime,functionbprime,B},
makeBesselTable[bessela,functiona,r+1,sqrtq/emu,qsign,q];
makeBesselTable[besselb,functionb,r+1,sqrtq*emu,qsign,q];
functionaprime=besselprime[bessela,functiona,qsign];
functionbprime=besselprime[besselb,functionb,qsign];
(1/(B[l,r,q,j,0])) mySum[Hold[modifieda*B[m,r,q,j,0] ((sqrtq*emu*functiona[m-l] functionbprime[m+l+p])-(sqrtq/emu functionb[m+l+p] functionaprime[m-l])-modifiedb*((sqrtq emu functiona[m+l+p] functionbprime[m-l])-(sqrtq/emu functionb[m-l] functionaprime[m+l+p])))],m,l,q]
]

(*Specific instances of the radial Mathieu functions and their derivatives*)
JE[r_,q_,\[Mu]_]:=KroneckerDelta[r] /; q==0
JO[r_,q_,\[Mu]_]:=0 /; q==0
JE[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RaHOCoE[r] RadialE[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]/;Re[q]>=0
JO[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RadialO[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]/;Re[q]>=0
JEPrime[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RaHOCoE[r] RadialEPrime[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]/;Re[q]>=0
JOPrime[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RadialOPrime[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]/;Re[q]>=0

IE[r_?EvenQ,q_,\[Mu]_]:=RaHOCoE[r] RadialE[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,Abs[r]/2,0,-1]
IO[r_?EvenQ,q_,\[Mu]_]:=RadialO[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,Abs[r]/2,0,-1]
IEPrime[r_?EvenQ,q_,\[Mu]_]:=RaHOCoE[r] RadialEPrime[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,Abs[r]/2,0,-1]
IOPrime[r_?EvenQ,q_,\[Mu]_]:=RadialOPrime[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,Abs[r]/2,0,-1]

YE[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RaHOCoE[r] RadialE[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]
YO[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RadialO[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]
YEPrime[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RaHOCoE[r] RadialEPrime[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]
YOPrime[r_?EvenQ,q_,\[Mu]_]:=RaNeg[r]RadialOPrime[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+Abs[r]/2),1,Abs[r]/2,0,1]

KE[r_?EvenQ,q_,\[Mu]_]:=RaPNCoE[r]RaHOCoE[r] RadialE[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,1,Abs[r]/2,0,-1]
KO[r_?EvenQ,q_,\[Mu]_]:=RaPNCoE[r] RadialO[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,1,Abs[r]/2,0,-1]
KEPrime[r_?EvenQ,q_,\[Mu]_]:=RaPNCoE[r]RaHOCoE[r]RadialEPrime[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,1,Abs[r]/2,0,-1]
KOPrime[r_?EvenQ,q_,\[Mu]_]:=RaPNCoE[r] RadialOPrime[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,1,Abs[r]/2,0,-1]

JE[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialE[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]/;Re[q]>=0
JO[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialO[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]/;Re[q]>=0
JEPrime[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialEPrime[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]/;Re[q]>=0
JOPrime[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialOPrime[Abs[r],q,\[Mu],BesselJ,BesselJ,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]/;Re[q]>=0

YE[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialE[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]
YO[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialO[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]
YEPrime[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialEPrime[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]
YOPrime[r_?OddQ,q_,\[Mu]_]:=RaNeg[r]RadialOPrime[Abs[r],q,\[Mu],BesselJ,BesselY,(-1)^(m+(Abs[r]-1)/2),1,(Abs[r]-1)/2,1,1]

IE[r_,q_,\[Mu]_]:=KroneckerDelta[r] /; q==0
IO[r_,q_,\[Mu]_]:=0 /; q==0
IE[r_?OddQ,q_,\[Mu]_]:=RadialE[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,(Abs[r]-1)/2,1,-1]
IO[r_?OddQ,q_,\[Mu]_]:=RadialO[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,(Abs[r]-1)/2,1,-1]
IEPrime[r_?OddQ,q_,\[Mu]_]:=RadialEPrime[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,(Abs[r]-1)/2,1,-1]
IOPrime[r_?OddQ,q_,\[Mu]_]:=RadialOPrime[Abs[r],-q,\[Mu],BesselI,BesselI,1,1,(Abs[r]-1)/2,1,-1]

KE[r_?OddQ,q_,\[Mu]_]:=RaPNCoER[r] RadialE[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,-1,(Abs[r]-1)/2,1,-1]
KO[r_?OddQ,q_,\[Mu]_]:=RaPNCoER[r] RadialO[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,-1,(Abs[r]-1)/2,1,-1]
KEPrime[r_?OddQ,q_,\[Mu]_]:=RaPNCoER[r] RadialEPrime[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,-1,(Abs[r]-1)/2,1,-1]
KOPrime[r_?OddQ,q_,\[Mu]_]:=RaPNCoER[r] RadialOPrime[Abs[r],-q,\[Mu],BesselI,BesselK,(-1)^m,-1,(Abs[r]-1)/2,1,-1]

JE[r_,q_,\[Mu]_]:=(I^r)IE[Abs[r],-q,\[Mu]]/;Re[q]<0
JO[r_,q_,\[Mu]_]:=(I^r)IO[Abs[r],-q,\[Mu]]/;Re[q]<0
JEPrime[r_,q_,\[Mu]_]:=(I^Abs[r])IEPrime[Abs[r],-q,\[Mu]]/;Re[q]<0
JOPrime[r_,q_,\[Mu]_]:=(I^Abs[r])IOPrime[Abs[r],-q,\[Mu]]/;Re[q]<0

(*special case for q=zero for radial Mathieu functions*)
makeRadialZeroFunc[name_]:=
(MakeExpression[RowBox[{ToString[name],"[r_,q_?isZero,\[Mu]_]",":=","0"}],StandardForm])

ReleaseHold[makeRadialZeroFunc/@shorthandFN[[5;;12;;1]]]


(*Series expansion of characteristic values*)

(*Building blocks for even characteristic values*)
\[CapitalAlpha][0,0,r_]:=1
\[CapitalAlpha][k_,0,r_]:=0
\[CapitalAlpha][k_,i_,r_?OddQ]:=If[Abs[i/2]>k,0,If[i==-2r,0,If[i==2k,((-1)^k)/((4^k)k!*Pochhammer[r+1,k]),If[i==-2k&&k<r,((-1)^k)/((4^k)k!*Pochhammer[1-r,k]),(1/(2i*(r+(i/2))))*(-\[CapitalAlpha][k-1,i+2,r]-\[CapitalAlpha][k-1,i-2,r]+If[(k+1-Abs[r+(i/2)])/2<1,0,Sum[\[Alpha][r,2j-1]\[CapitalAlpha][k-2j+1,-i-2r,r],{j,1,(k+1-Abs[r+(i/2)])/2}]]+If[((k-Abs[i/2])/2)<1,0,Sum[\[Alpha][r,2j]*\[CapitalAlpha][k-2j,i,r],{j,1,((k-Abs[i/2])/2)}]])]]]]
\[CapitalAlpha][k_,i_,r_?EvenQ]:=If[Abs[i/2]>k,0,If[i==-2r,0,If[i==2k,((-1)^k)/((4^k)k!*Pochhammer[r+1,k]),If[i==-2k&&k<r,((-1)^k)/((4^k)k!*Pochhammer[1-r,k]),(1/(2i(r+(i/2))))(-\[CapitalAlpha][k-1,i+2,r]-\[CapitalAlpha][k-1,i-2,r]+If[((k-Abs[i/2])/2)<1,0,Sum[\[Alpha][r,2j]\[CapitalAlpha][k-2j,i,r],{j,1,((k-Abs[i/2])/2)}]])]]]]

(*Even characteristic values*)
a[ r_,q_,n_]:=a[r,q,n]=r^2+Sum[\[Alpha][Abs[r],k]q^k,{k,1,n}]
\[Alpha][r_?OddQ,k_?OddQ]:=\[Alpha][r,k]=\[CapitalAlpha][k-1,2-2r,r]+\[CapitalAlpha][k-1,-2-2r,r]
\[Alpha][r_?EvenQ,k_?OddQ]:=0
\[Alpha][r_?EvenNotZero,k_?EvenQ]:=\[Alpha][r,k]=\[CapitalAlpha][k-1,2,r]+\[CapitalAlpha][k-1,-2,r]+\[CapitalAlpha][k-1,2-2r,r]+\[CapitalAlpha][k-1,-2-2r,r]
\[Alpha][r_?OddZero,k_?EvenQ]:=\[Alpha][r,k]=\[CapitalAlpha][k-1,2,r]+\[CapitalAlpha][k-1,-2,r]

(*Building blocks for odd characteristic values*)
\[CapitalBeta][0,0,r_]:=1
\[CapitalBeta][k_,0,r_]:=0
\[CapitalBeta][k_,i_,r_?OddQ]:=If[Abs[i/2]>k,0,If[i==-2r,0,If[i==2k,((-1)^k)/((4^k)k!*Pochhammer[r+1,k]),If[i==-2k&&k<r,((-1)^k)/((4^k)k!*Pochhammer[1-r,k]),(1/(2i*(r+(i/2))))*(-\[CapitalBeta][k-1,i+2,r]-\[CapitalBeta][k-1,i-2,r]-If[(k+1-Abs[r+(i/2)])/2<1,0,Sum[(\[Beta][r,2j-1]\[CapitalBeta][k-2j+1,-i-2r,r]),{j,1,(k+1-Abs[r+(i/2)])/2}]]+If[((k-Abs[i/2])/2)<1,0,Sum[\[Beta][r,2j]\[CapitalBeta][k-2j,i,r],{j,1,(k-Abs[i/2])/2}]])]]]]
\[CapitalBeta][k_,i_,r_?EvenQ]:=If[Abs[i/2]>k,0,If[i==-2r,0,If[i==2k,((-1)^k)/((4^k)k!*Pochhammer[r+1,k]),If[i==-2k&&k<r,((-1)^k)/((4^k)k!*Pochhammer[1-r,k]),(1/(2i*(r+(i/2))))*(-\[CapitalBeta][k-1,i+2,r]-\[CapitalBeta][k-1,i-2,r]+If[((k-Abs[i/2])/2)<1,0,Sum[\[Beta][r,2j]\[CapitalBeta][k-2j,i,r],{j,1,((k-Abs[i/2])/2)}]])]]]]

(*Odd characteristic values*)
b[ r_,q_,n_]:=b[r,q,n]=r^2+Sum[\[Beta][Abs[r],k]q^k,{k,1,n}]
b[r_?isZero,q_,n_]:=Message[MathieuCharacteristicB::zord]
\[Beta][r_?OddQ,k_?OddQ]:=\[Beta][r,k]=-\[CapitalBeta][k-1,2-2r,r]+-\[CapitalBeta][k-1,-2-2r,r]
\[Beta][r_?EvenQ,k_?OddQ]:=0
\[Beta][r_?EvenQ,k_?EvenQ]:=\[Beta][r,k]=\[CapitalBeta][k-1,2,r]+\[CapitalBeta][k-1,-2,r]-\[CapitalBeta][k-1,2-2r,r]-\[CapitalBeta][k-1,-2-2r,r]
\[Beta][r_?OddQ,k_?EvenQ]:=\[Beta][r,k]=\[CapitalBeta][k-1,2,r]+\[CapitalBeta][k-1,-2,r]

(*Summations for Series function*)
ce[r_?OddQ,q_,theta_,i_]:=(1/Sum[(A[m,Abs[r],q,i+1,1])^2,{m,0,Floor[Abs[r]/2]+i}]^(1/2))Sum[A[m,Abs[r],q,i+1,1]Cos[(2m+1)theta],{m,0,Floor[Abs[r]/2]+i}]
ce[r_?EvenQ,q_,theta_,i_]:=(1/(1+Sum[A[m,Abs[r],q,i+1,1]^2,{m,0,(Abs[r]/2)+i}])^(1/2))Sum[A[m,Abs[r],q,i+1,1]Cos[(2m)theta],{m,0,(Abs[r]/2)+i}]
se[r_?OddQ,q_,theta_,i_]:=If[Negative[r],-1,1](1/Sum[(B[m,Abs[r],q,i+1,1])^2,{m,0,Floor[Abs[r]/2]+i}]^(1/2))Sum[B[m,Abs[r],q,i+1,1]Sin[(2m+1)theta],{m,0,Floor[Abs[r]/2]+i}]
se[r_?EvenQ,q_,theta_,i_]:=If[Negative[r],-1,1](1/(Sum[B[m,Abs[r],q,i+1,1]^2,{m,0,(Abs[r]/2)+i}])^(1/2))Sum[B[m,Abs[r],q,i+1,1]Sin[(2m)theta],{m,0,(Abs[r]/2)+i}]  
je[r_?EvenQ,q_,mu_,i_]:=RaNeg[r]RaHOCoE[r](1/(A[Abs[r]/2,Abs[r],q,i,1])) Sum[(-1)^(m + Abs[r]/2) A[m,Abs[r],q,i,1](BesselJ[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselJ[m+Abs[r]/2,Sqrt[q]Exp[mu]] +  BesselJ[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselJ[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]/2-Floor[i/2]>0,Floor[(Abs[r]-i+2)/2],0],Floor[(Abs[r]/2+i)/2]}]
je[r_?OddQ,q_,mu_,i_]:=RaNeg[r](1/(A[(Abs[r]-1)/2,Abs[r],q,i,1])) Sum[(-1)^(m + (Abs[r]-1)/2) A[m,Abs[r],q,i,1](BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] +  BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]/2-Floor[(i+1)/2]>0,Floor[(Abs[r]-i+2)/2],0],Floor[(Abs[r]/2+i)/2]-1}]
ie[r_?EvenQ, q_, mu_,i_]:=RaHOCoE[r](1/A[Abs[r]/2,Abs[r],-q,i,1]) Sum[A[m,Abs[r],-q,i,1](BesselI[m-Abs[r]/2,Sqrt[q]Exp[-mu]]BesselI[m+Abs[r]/2,Sqrt[q]Exp[mu]]+BesselI[m-Abs[r]/2,Sqrt[q]Exp[mu]]BesselI[m+Abs[r]/2,Sqrt[q]Exp[-mu]]),{m,If[Abs[r]/2-Floor[i/2]>0,Floor[(Abs[r]-i+2)/2],0], Floor[(Abs[r]/2+1)/2]}]
ie[r_?OddQ, q_, mu_,i_]:=(1/A[(Abs[r]-1)/2,Abs[r],-q,i,1]) Sum[ A[m,Abs[r],-q,i,1](BesselI[m-(Abs[r]-1)/2,Sqrt[q] Exp[-mu]]BesselI[m+(Abs[r]-1)/2+1,Sqrt[q] E^mu]+BesselI[m-(Abs[r]-1)/2,Sqrt[q] E^mu]BesselI[m+(Abs[r]-1)/2+1,Sqrt[q] E^-mu]),{m,If[Abs[r]/2-Floor[(i+1)/2]>0,Floor[(Abs[r]-i+2)/2],0],Floor[(Abs[r]/2+i-1)/2]}]
jo[r_?EvenQ,q_,mu_,i_]:=RaNeg[r](1/(B[Abs[r]/2,Abs[r],q,i,1])) Sum[(-1)^(m + Abs[r]/2) B[m,Abs[r],q,i,1](BesselJ[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselJ[m+Abs[r]/2,Sqrt[q]Exp[mu]] -  BesselJ[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselJ[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]-2i>0,Abs[r]/2-i+1,1],Floor[(Abs[r]+i-1)/2]}]
jo[r_?OddQ,q_,mu_,i_]:=RaNeg[r](1/(B[(Abs[r]-1)/2,Abs[r],q,i,1])) Sum[(-1)^(m + (Abs[r]-1)/2) B[m,Abs[r],q,i,1](BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] -  BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m, If[Abs[r]-2i>0,Floor[(Abs[r]-2i+1)/2],0],Floor[(Abs[r]+i-1)/2]}]
io[r_?EvenQ, q_, mu_,i_]:=(1/(B[Abs[r]/2,Abs[r],-q,i,1])) Sum[ B[m,Abs[r],-q,i,1](BesselI[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselI[m+Abs[r]/2,Sqrt[q]Exp[mu]] -  BesselI[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselI[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]/2-Floor[(i+1)/2]>0,Floor[(Abs[r]-i+2)/2]+1,1],Floor[(Abs[r]/2+i)/2]}]
io[r_?OddQ, q_, mu_,i_]:=(1/(B[(Abs[r]-1)/2,Abs[r],-q,i,1])) Sum[ B[m,Abs[r],-q,i,1](BesselI[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselI[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] -  BesselI[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselI[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]/2-Floor[(i+1)/2]>0,Floor[(Abs[r]-i+2)/2],0],Floor[(Abs[r]/2+i)/2]}]
ye[r_?EvenQ, q_, mu_,i_] :=RaNeg[r]RaHOCoE[r] (1/(A[Abs[r]/2,Abs[r],q,i,1])) Sum[(-1)^(m + Abs[r]/2) A[m,Abs[r],q,i,1](BesselJ[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselY[m+Abs[r]/2,Sqrt[q]Exp[mu]] +  BesselJ[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselY[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+2, Floor[(Abs[r]-4i)/2],0],Abs[r]/2+i+1}]
ye[r_?OddQ, q_, mu_,i_] :=RaNeg[r](1/(A[(Abs[r]-1)/2,Abs[r],q,i,1])) Sum[(-1)^(m + (Abs[r]-1)/2) A[m,Abs[r],q,i,1](BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselY[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] +  BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselY[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m, 0,(Abs[r]+2i+1)/2}]
ke[r_?EvenQ, q_, mu_,i_] :=RaPNCoE[r]RaHOCoE[r](1/(A[Abs[r]/2,Abs[r],-q,i,1])) Sum[(-1)^m A[m,Abs[r],-q,i,1](BesselI[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselK[m+Abs[r]/2,Sqrt[q]Exp[mu]] +  BesselI[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselK[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+2, Floor[(Abs[r]-4i)/2],0],If[Abs[r]>1,2+i+Floor[(Abs[r]-2)/2],i+3]}]
ke[r_?OddQ, q_, mu_,i_] :=RaPNCoER[r] (1/(A[(Abs[r]-1)/2,Abs[r],-q,i,1])) Sum[ (-1)^m A[m,Abs[r],-q,i,1](BesselI[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselK[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] - BesselI[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselK[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+2, Floor[(Abs[r]-4i)/2],0],If[Abs[r]>1,2+i+Floor[(Abs[r]-2)/2],i+3]}]
yo[r_?EvenQ, q_, mu_,i_] :=RaNeg[r](1/(B[Abs[r]/2,Abs[r],q,i,1])) Sum[(-1)^(m + Abs[r]/2) B[m,Abs[r],q,i,1](BesselJ[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselY[m+Abs[r]/2,Sqrt[q]Exp[mu]] +  BesselJ[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselY[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+2, Floor[(Abs[r]-4i)/2],0],Abs[r]/2+i+1}]
yo[r_?OddQ, q_, mu_,i_] :=RaNeg[r](1/(B[(Abs[r]-1)/2,Abs[r],q,i,1])) Sum[(-1)^(m + (Abs[r]-1)/2) B[m,Abs[r],q,i,1](BesselJ[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselY[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] +  BesselJ[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselY[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m, 0,(Abs[r]+2i+1)/2}]
ko[r_?EvenQ, q_, mu_,i_] :=RaPNCoE[r] (1/(B[Abs[r]/2,Abs[r],-q,i,1])) Sum[(-1)^m B[m,Abs[r],-q,i,1](BesselI[m-Abs[r]/2,Sqrt[q]Exp[-mu]] BesselK[m+Abs[r]/2,Sqrt[q]Exp[mu]] +  BesselI[m+Abs[r]/2,Sqrt[q]Exp[-mu]] BesselK[m-Abs[r]/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+4, Floor[(Abs[r]-4i)/2],0],If[Abs[r]>1,2+i+Floor[(Abs[r]-2)/2],i+3]}]
ko[r_?OddQ, q_, mu_,i_] := RaPNCoER[r] (1/(B[(Abs[r]-1)/2,Abs[r],-q,i,1])) Sum[(-1)^m B[m,Abs[r],-q,i,1](BesselI[m-(Abs[r]-1)/2,Sqrt[q]Exp[-mu]] BesselK[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[mu]] - BesselI[m+(Abs[r]-1)/2+1,Sqrt[q]Exp[-mu]] BesselK[m-(Abs[r]-1)/2,Sqrt[q]Exp[mu]]), {m,If[Abs[r]>=4i+4, Floor[(Abs[r]-4i)/2],0],If[Abs[r]>1,2+i+Floor[(Abs[r]-2)/2],i+3]}]


Unprotect[Series];

makeSeries[func1_,func2_] := (Series[func1[r_,z_,theta_],{q_,0,n_}]:=Simplify[Series[func2[r,z,theta,n],{q,0,n}],{theta\[Element]Reals}]/;!FreeQ[z,q])
makeSeriesDeriv[func1_,func2_] := (Series[func1[r_,z_,theta_],{q_,0,n_}]:=D[Series[func2[r,z,theta],{q,0,n}],theta])

makeSeries[MathieuCE,ce]
makeSeries[MathieuSE,se]
makeSeries[MathieuJE,je]
makeSeries[MathieuJO,jo]
makeSeries[MathieuYE,ye]
makeSeries[MathieuYO,yo]
makeSeries[MathieuIE,ie]
makeSeries[MathieuIO,io]
makeSeries[MathieuKE,ke]
makeSeries[MathieuKO,ko]
makeSeriesDeriv[MathieuCEPrime, MathieuCE]
makeSeriesDeriv[MathieuSEPrime, MathieuSE]
makeSeriesDeriv[MathieuJEPrime, MathieuJE]
makeSeriesDeriv[MathieuJOPrime, MathieuJO]
makeSeriesDeriv[MathieuYEPrime, MathieuYE]
makeSeriesDeriv[MathieuYOPrime, MathieuYO]
makeSeriesDeriv[MathieuIOPrime, MathieuIO]
makeSeriesDeriv[MathieuIEPrime, MathieuIE]
makeSeriesDeriv[MathieuKEPrime, MathieuKE]
makeSeriesDeriv[MathieuKOPrime, MathieuKO]
Series[MathieuCharacteristicA[r_,z_],{q_,0,n_}]:=Series[a[r,z,n],{q,0,n}]/;!FreeQ[z,q]
Series[MathieuCharacteristicB[r_,z_],{q_,0,n_}]:=Series[b[r,z,n],{q,0,n}]/;!FreeQ[z,q]
Protect[Series];


clearValues[f_,n_:1]:=DownValues[f]=Take[DownValues[f],-n]

MathieuClearCache[]:=(clearValues[Mathieu`Private`Ae,2];clearValues[Mathieu`Private`Ao,2];clearValues[Mathieu`Private`Be,3];clearValues[Mathieu`Private`Bo,2];
clearValues[Mathieu`Private`Te,4];clearValues[Mathieu`Private`To,3];
clearValues[Mathieu`Private`G];clearValues[Mathieu`Private`Qe,3];clearValues[Mathieu`Private`Qo,3];
clearValues[Mathieu`Private`JFOdd];clearValues[Mathieu`Private`JFEven];
Null)

End[]

SetAttributes[{MathieuCE,MathieuSE,MathieuFE,MathieuFO,MathieuJE,MathieuJO,MathieuYE,MathieuYO,MathieuIE,MathieuIO,MathieuKE,MathieuKO,MathieuCEPrime,MathieuSEPrime,MathieuFEPrime,MathieuFOPrime,MathieuJEPrime,MathieuJOPrime,MathieuYEPrime,MathieuYOPrime,MathieuIEPrime,MathieuIOPrime,MathieuKEPrime,MathieuKOPrime,MathieuAngularEvenWronskian,MathieuAngularOddWronskian,MathieuRadialEvenWronskian,MathieuRadialOddWronskian,MathieuModifiedEvenWronskian,MathieuModifiedOddWronskian},{Listable,Protected,ReadProtected,NumericFunction}];

EndPackage[]
