(*******************************************************************

This package is based on previous work of Falloon (see P. E. Falloon, P. C.
Abbott, and J. B.  Wang, J. Phys. A36 (2003) 5477), which was released
without license restrictions or copyright.  Modifications are copyright
&copy Noah Graham, 2005-2008, all rights reserved.  They may be used and
modified for noncommercial research purposes, provided that citation to 

T. Emig, N. Graham, R. L. Jaffe and M. Kardar,
"Casimir manipulations: The Orientation Dependence of
Fluctuation-Induced Forces," arXiv:0811.1597,
Phys. Rev. D77 (2008) 025005.

is included in all publications or other products in which the modified
program or any programs derived from it was used.  Commercial use or
inclusion in a commercial product is prohibited.  This program and all
accompanying materials are provided "AS IS," WITHOUT WARRANTY OF ANY KIND,
EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.

***********************************************************************)

(* bug workaround for series at Infinity*)
Off[SeriesData::sdatv]

Unprotect["Spheroidal*"]
ClearAll["Spheroidal*"]

BeginPackage["Spheroidal`"]; 

Unprotect["Spheroidal`*"]; 

ClearAll["Spheroidal`*", "Spheroidal`Private`*"]; 

SpheroidalEigenvalue::usage = "SpheroidalEigenvalue[\[Nu], \[Mu], \[Gamma]] is \
Meixner's spheroidal eigenvalue with degree \[Nu], order \[Mu], and spheroidicity \
\[Gamma]."; 
SpheroidalPS::usage = "SpheroidalPS[\[Nu], \[Mu], type:1, z, \[Gamma]] is Meixner's \
angular spheroidal function of the first kind, with possible types 1 or 2."; 
SpheroidalPSPrime::usage = "SpheroidalPSPrime[\[Nu], \[Mu], type:1, z, \[Gamma]] is \
the derivative of Meixner's angular spheroidal function of the first kind, \
with possible types 1 or 2."; 
SpheroidalQS::usage = "SpheroidalQS[\[Nu], \[Mu], type:1, z, \[Gamma]] is Meixner's \
angular spheroidal function of the second kind, with possible types 1 or 2."; 
SpheroidalQSPrime::usage = "SpheroidalQSPrime[\[Nu], \[Mu], type:1, z, \[Gamma]] is \
the derivative of Meixner's angular spheroidal function of the second kind, \
with possible types 1 or 2."; 
SpheroidalS1::usage = "SpheroidalS1[\[Nu], \[Mu], z, \[Gamma]] is Meixner's radial \
spheroidal function of the first kind."; 
SpheroidalS1Prime::usage = "SpheroidalS1Prime[\[Nu], \[Mu], z, \[Gamma]] is the \
derivative of Meixner's radial spheroidal function of the first kind."; 
SpheroidalS2::usage = "SpheroidalS2[\[Nu], \[Mu], z, \[Gamma]] is Meixner's radial \
spheroidal function of the second kind."; 
SpheroidalS2Prime::usage = "SpheroidalS2Prime[\[Nu], \[Mu], z, \[Gamma]] is the \
derivative of Meixner's radial spheroidal function of the second kind."; 
SpheroidalS3::usage = "SpheroidalS3[\[Nu], \[Mu], z, \[Gamma]] is Meixner's radial \
spheroidal function of the third kind."; 
SpheroidalS3Prime::usage = "SpheroidalS3Prime[\[Nu], \[Mu], z, \[Gamma]] is the \
derivative of Meixner's radial spheroidal function of the third kind."; 
SpheroidalRadialFactor::usage = "SpheroidalRadialFactor[\[Nu], \[Mu], \[Gamma]] is \
Meixner's radial spheroidal normalization factor."; 
SpheroidalJoiningFactor::usage = "SpheroidalJoiningFactor[\[Nu], \[Mu], \[Gamma]] is \
Meixner's spheroidal joining factor.";
SpheroidalClearCache::usage = "Clears internal caches to save memory."


Begin["`Private`"]; 

SetAttributes[{friendlyOff, friendlyOn, messageStatus}, HoldAll]; 
messageStatus[MessageName[func_, mesg_String]] := 
  If[Head[MessageName[func, mesg] /. Messages[func]] =!= $Off, $On, $Off, 
   $Off]
friendlyOff[mesg_MessageName] := 
  (mesgWasOn[Hold[mesg]] = messageStatus[mesg] =!= $Off; Off[mesg]; )
friendlyOn[mesg_MessageName] := If[mesgWasOn[Hold[mesg]], 
    On[mesg]; mesgWasOn[Hold[mesg]] = False; ]; 

friendlyOff[General::spell]; 
friendlyOff[General::spell1]; 

precRound[z_] := If[z==MachinePrecision, MachinePrecision, Max[MachinePrecision, 10 Ceiling[Round[z]/10]]]

numQ[z__] := And @@ NumericQ /@ {z}
intQ[n__] := numQ[n] && Round[Re[{n}]] == {n}
natQ[n__] := intQ[n] && And @@ NonNegative[{n}]
eltQ[z_, list_] := numQ[z] && Or @@ (z == #1 & ) /@ list
zerQ[z__] := numQ[z] && And @@ (#1 == 0 & ) /@ {z}
memQ[{lower_, upper_}, z_] := IntervalMemberQ[Interval[{lower, upper}], 
   Chop[z]]

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

closestElt[list_List, num_] := Module[{dist, minpos}, 
   dist = Abs[list - num]; minpos = 
     First[Flatten[Position[dist, Min[dist]]]]; list[[minpos]]]

zerorule = {x_ /; x == 0. -> 0}

legP[nu_, 0, a_, z_] := LegendreP[nu, z]
legP[nu_, mu_, 1, z_] := KroneckerDelta[mu] /; intQ[nu, mu] && z==1
legP[nu_, mu_, 2, z_] := KroneckerDelta[mu] /; intQ[nu, mu] && z==1

legP[nu_, mu_, 1, z_] := ((1 + z)^(mu/2)/(1 - z)^(mu/2))*
   Hypergeometric2F1Regularized[-nu, nu + 1, 1 - mu, (1 - z)/2]

legP[nu_, mu_, 2, z_] := ((z + 1)^(mu/2)/(z - 1)^(mu/2))*
   Hypergeometric2F1Regularized[-nu, nu + 1, 1 - mu, (1 - z)/2]

legQ[nu_, mu_, 2, z_] /; numQ[z] &&  !memQ[{-1, 1}, z] := 
  2^(-nu - 1)*Exp[I*mu*Pi]*Sqrt[Pi]*Gamma[mu + nu + 1]*z^(-mu - nu - 1)*
   (z + 1)^(mu/2)*(z - 1)^(mu/2)*Hypergeometric2F1Regularized[
    (mu + nu + 1)/2, (mu + nu + 2)/2, nu + 3/2, 1/z^2]

legQ[nu_, mu_, 1, z_] /; numQ[z] &&  !memQ[{-1, 1}, z] &&  !intQ[mu] := 
  ((1 - z)^(mu/2)/(z - 1)^(mu/2))*(legQ[nu, mu, 2, z]/E^(I*mu*Pi) + 
    (Pi/2)*Csc[mu*Pi]*(((z - 1)^mu/(1 - z)^mu)*Cos[mu*Pi] - 1)*
     legP[nu, mu, 2, z])

legQ[nu_, m_, 1, z_] /; numQ[z] &&  !memQ[{-1, 1}, z] && intQ[m] := 
  (-1)^m*((1 - z)^(m/2)/(z - 1)^(m/2))*(legQ[nu, m, 2, z] + 
    (Pi/2)*(Sqrt[z - 1]/Sqrt[1 - z])*legP[nu, m, 2, z])

legQ[nu_, mu_, 1, z_] /; memQ[{-1, 1}, z] := LegendreQ[nu, mu, 2, z]

legQ[nu_, mu_, 2, z_] /; memQ[{-1, 1}, z] &&  !intQ[mu] := 
  E^(I*mu*Pi)*((z - 1)^(mu/2)/(1 - z)^(mu/2))*(legQ[nu, mu, 1, z] - 
    (Pi/2)*Csc[mu*Pi]*(Cos[mu*Pi] - (1 - z)^mu/(z - 1)^mu)*legP[nu, mu, 1, z])

legQ[nu_, m_, 2, z_] /; memQ[{-1, 1}, z] && intQ[m] := 
  (-1)^m*((z - 1)^(m/2)/(1 - z)^(m/2))*legQ[nu, m, 1, z] - 
   (1/2)*(Log[z - 1] - Log[1 - z])*legP[nu, m, 2, z]

legPPrime[nu_, mu_, a_, z_] := (1/(z^2 - 1))*
   ((-(mu + nu))*legP[nu - 1, mu, a, z] + nu*z*legP[nu, mu, a, z])

legQPrime[nu_, mu_, a_, z_] /;  !nu == -mu := 
  (1/(z^2 - 1))*((-(mu + nu))*legQ[nu - 1, mu, a, z] + 
    nu*z*legQ[nu, mu, a, z])

legQPrime[nu_, mu_, 1, z_] /; nu == -mu := 
  D[LegendreQ[Rationalize[nu], Rationalize[mu], zz], zz] /. zz -> z

legQPrime[nu_, mu_, 2, z_] /; nu == -mu &&  !intQ[mu] := 
  Exp[I*mu*Pi]*((z - 1)^(mu/2)/(1 - z)^(mu/2))*(legQPrime[nu, mu, 1, z] - 
    (Pi/2)*Csc[mu*Pi]*(Cos[mu*Pi] - (z - 1)^mu/(1 - z)^mu)*
     legPPrime[nu, mu, 1, z])

legQPrime[n_, m_, 2, z_] /; intQ[m, n] && n == -m := 
  (-1)^m*((z - 1)^(m/2)/(1 - z)^(m/2))*(legQPrime[n, m, 1, z] - 
    (Pi/2)*(Sqrt[z - 1]/Sqrt[1 - z])*legPPrime[n, m, 1, z])

besJ[nu_, z_] := Sqrt[Pi/2]*(BesselJ[nu + 1/2, z]/Sqrt[z])
besY[nu_, z_] := Sqrt[Pi/2]*(BesselY[nu + 1/2, z]/Sqrt[z])
besK[nu_, z_] := Sqrt[2/Pi]*(BesselK[nu + 1/2, z/I]/((I^(nu+2))*Sqrt[z/I]))

$precInc = 50;
$precIncBig = 100;
$precIncInit = 20;
$precIncInitSmall = 10;
$seriesInc = 40;

aLS[nu_, mu_] = ((mu - nu - 2)*(mu - nu - 1))/((2*nu + 1)*(2*nu + 3)); 
bLS[nu_, mu_] = (2*nu*(nu + 1) - 2*mu*mu - 1)/((2*nu - 1)*(2*nu + 3)); 
cLS[nu_, mu_] = ((mu + nu - 1)*(mu + nu))/((2*nu - 1)*(2*nu + 1)); 


legTable[nu_, mu_, type_, kmin_, kmax_, z_, leg_] := 
  Module[{prec, nu1, mu1, z1, list}, prec = Precision[z] + $precIncBig;
    While[True, {nu1, mu1} = SetPrecision[{nu, mu}, 
        Max[prec, Precision[{nu, mu}]]]; z1 = SetPrecision[z, prec]; 
      list = {leg[nu1, mu1, type, z1], leg[nu1 + 2, mu1, type, z1]}; 
      Do[AppendTo[list, If[aLS[nu1 + 2*k, mu1] == 0, leg[nu1 + 2*k + 2, 
          mu1, type, z1], -((bLS[nu1 + 2*k, mu1] - z1^2)*list[[-1]] + 
            cLS[nu1 + 2*k, mu1]*list[[-2]])/aLS[nu1 + 2*k, mu1]]], 
       {k, 1, kmax - 1}]; Do[PrependTo[list, If[cLS[nu1 + 2*k, mu1] == 0, 
         leg[nu1 + 2*k - 2, mu1, type, z1], 
         -((bLS[nu1 + 2*k, mu1] - z1^2)*list[[1]] + aLS[nu1 + 2*k, mu1]*
             list[[2]])/cLS[nu1 + 2*k, mu1]]], {k, 0, kmin + 1, -1}]; 
      list = list /. zerorule; If[Precision[list] >= Precision[z], 
       Break[], prec += $precInc]]; N[list, Precision[z]]]

legPTable[nu_, mu_, type_, kmin_, kmax_, z_] := legTable[nu, mu, type, kmin, kmax, z, legP]
legQTable[nu_, mu_, type_, kmin_, kmax_, z_] := legTable[nu, mu, type, kmin, kmax, z, legQ]

mypoch[a_,n_]:= Module[{j, ret=1.0},For[j=0, j<n, j++, ret *= (j+a)]; ret]

legPRecursionCompiledReal = Compile[{
    {nu,_Integer},{mu,_Integer},{m,_Integer},{z1,_Real},{startindex,_Integer},{startp0,_Real},{startp1,_Real},{kmin,_Integer},{kmax,_Integer}},
    Module[{index=startindex,listindex=1,p0=startp0,p1=startp1,list=Table[0.,{kmax-kmin+1}]},
    While[True,
         list[[listindex++]] = If[m<0, If[Mod[mu,2]==1,-1,1]*p0/mypoch[1.0*(index-mu),2.0*mu], p0];
         If[index < nu+2*kmax,
            p0 = (z1*(2*index+1)*p1 - (index+mu)*p0)/(index+1-mu);
            p1 = (z1*(2*index+3)*p0 - (index+mu+1)*p1)/(index+2-mu); index = index + 2, Break[]]]; list]]

legPRecursion[nu_Integer,mu_Integer,m_Integer,z1_Real,index_Integer,p0_Real,p1_Real,kmin_Integer,kmax_Integer] := legPRecursionCompiledReal[nu,mu,m,z1,index,p0,p1,kmin,kmax] /; Precision[z1]==MachinePrecision


legPRecursionCompiledComplex = Compile[{
    {nu,_Integer},{mu,_Integer},{m,_Integer},{z1,_Complex},{startindex,_Integer},{startp0,_Complex},{startp1,_Complex},{kmin,_Integer},{kmax,_Integer}},
    Module[{index=startindex,listindex=1,p0=startp0,p1=startp1,list=Table[0.+0. I,{kmax-kmin+1}]},
    While[True,
         list[[listindex++]] = If[m<0, If[Mod[mu,2]==1,-1,1]*p0/mypoch[1.0*(index-mu),2.0*mu], p0];
         If[index < nu+2*kmax,
            p0 = (z1*(2*index+1)*p1 - (index+mu)*p0)/(index+1-mu);
            p1 = (z1*(2*index+3)*p0 - (index+mu+1)*p1)/(index+2-mu); index = index + 2, Break[]]]; list]]


legPRecursion[nu_Integer,mu_Integer,m_Integer,z1_Complex,index_Integer,p0_Complex,p1_Complex,kmin_Integer,kmax_Integer] := legPRecursionCompiledComplex[nu,mu,m,z1,index,p0,p1,kmin,kmax] /; Precision[z1]==MachinePrecision

legPRecursion[nu_,mu_,m_,z1_,startindex_,startp0_,startp1_,kmin_,kmax_] :=
    Module[{index=startindex,listindex=1,p0=startp0,p1=startp1,list=Table[0,{kmax-kmin+1}]},
    While[True,
         list[[listindex++]] = If[m<0, If[Mod[mu,2]==1,-1,1]*p0/Pochhammer[index-mu,2*mu], p0];
         If[index < nu+2*kmax,
            p0 = (z1*(2*index+1)*p1 - (index+mu)*p0)/(index+1-mu);
            p1 = (z1*(2*index+3)*p0 - (index+mu+1)*p1)/(index+2-mu); index = index + 2, Break[]]]; list]

legPTableNoCheck[nu_, m_, 1, kmin_, kmax_, z_] :=
     Module[{z1 = z, mu=Abs[m], list, p0, p1, p2, index},
     While[True,
       p0 = If[z1*z1 == 1, KroneckerDelta[mu], (-1)^mu*((2*mu-1)!!)*(1 - z1*z1)^(mu/2)];
       p1 = z1*(2*mu+1)*p0;
       index = mu+1;

       While[index <= nu+2*kmin,
            p2 = (z1*(2*index+1)*p1 - (index+mu)*p0)/(index+1-mu); p0 = p1; p1 = p2; index++];

       list = legPRecursion[nu,mu,m,z1,index,p0,p1,kmin,kmax] /. zerorule;

       If[Precision[list] >= Precision[z], Break[], z1 = SetPrecision[z1, Precision[z1] + $precInc]]];
       N[list, Precision[z]]]

legPTableNoCheck[nu_, mu_?EvenQ, 2, kmin_, kmax_, z_] :=
  legPTableNoCheck[nu, mu, 1, kmin, kmax, z]*(-1)^(mu/2)

legPTableNoCheck[nu_, mu_, 2, kmin_, kmax_, z_] := legPTable[nu, mu, 2, kmin, kmax, z] /; z==1

legPTableNoCheck[nu_, mu_, 2, kmin_, kmax_, z_] :=
  legPTableNoCheck[nu, mu, 1, kmin, kmax, z]*(1-z)^(mu/2)/(z-1)^(mu/2)

legPTableNoCheck[nu_, mu_, type_, kmin_, kmax_, z_] := legPTable[nu, mu, type, kmin, kmax, z]

  
legPrimeTable[nu_, mu_, type_, kmin_, kmax_, z_, legTable_] := 
  Module[{nu1, mu1}, {nu1, mu1} = SetPrecision[{nu, mu}, Max[Precision[{nu, mu}], Precision[z]]]; 
    ((-Range[mu1 + nu1 + 2*kmin, mu1 + nu1 + 2*kmax, 2])*
       legTable[nu1 - 1, mu1, type, kmin, kmax, z] + 
      z*Range[nu1 + 2*kmin, nu1 + 2*kmax, 2]*legTable[nu1, mu1, type, kmin, 
        kmax, z])/(z^2 - 1)]

legPPrimeTable[nu_, mu_, type_, kmin_, kmax_, z_] := legPrimeTable[nu, mu, type, kmin, kmax, z, legPTable]

legQPrimeTable[nu_, mu_, type_, kmin_, kmax_, z_] /; intQ[mu + nu] := 
  Module[{nu1, mu1}, {nu1, mu1} = SetPrecision[{nu, mu}, Max[Precision[{nu, mu}], Precision[z]]]; 
    Prepend[((-Range[mu1 + nu1 + 2*kmin + 2, mu1 + nu1 + 2*kmax, 2])*
        legQTable[nu1 + 1, mu1, type, kmin, kmax - 1, z] + 
       z*Range[nu1 + 2*kmin + 2, nu1 + 2*kmax, 2]*legQTable[nu1 + 2, mu1, 
         type, kmin, kmax - 1, z])/(z^2 - 1), legQPrime[nu1 + 2*kmin, mu1, 
      type, z]]]

legQPrimeTable[nu_, mu_, type_, kmin_, kmax_, z_] /;  !intQ[mu + nu] :=
   legPrimeTable[nu, mu, type, kmin, kmax, z, legQTable]

aBS[nu_] = 1/((2*nu + 1)*(2*nu + 3)); 
bBS[nu_] = 2/((2*nu - 1)*(2*nu + 3));
cBS[nu_] = 1/((2*nu - 1)*(2*nu + 1));

besTable[nu_, kmin_, kmax_, z_, bes_] := Table[bes[nu+2*k, z], {k,kmin,kmax}]

besJTable[nu_, kmin_, kmax_, z_] := besTable[nu, kmin, kmax, z, besJ]
besYTable[nu_, kmin_, kmax_, z_] := besTable[nu, kmin, kmax, z, besY]
besKTable[nu_, kmin_, kmax_, z_] := besTable[nu, kmin, kmax, z, besK]
   
besPrimeTable[nu_, kmin_, kmax_, z_, bes_] := Table[Derivative[0,1][bes][nu+k, z], {k,2*kmin,2*kmax,2}]

besJPrimeTable[nu_, kmin_, kmax_, z_] := besPrimeTable[nu, kmin, kmax, z, besJ]
besYPrimeTable[nu_, kmin_, kmax_, z_] := besPrimeTable[nu, kmin, kmax, z, besY]
besKPrimeTable[nu_, kmin_, kmax_, z_] := besPrimeTable[nu, kmin, kmax, z, besK]
   
recA[nu_, mu_, k_, gsq_] = (-gsq)*(((nu - mu + 2*k - 1)*(nu - mu + 2*k))/
     ((2*nu + 4*k - 3)*(2*nu + 4*k - 1))); 
recB[nu_, mu_, k_, gsq_] = (nu + 2*k)*(nu + 2*k + 1) - 
    2*gsq*(((nu + 2*k)*(nu + 2*k + 1) + mu^2 - 1)/((2*nu + 4*k - 1)*
       (2*nu + 4*k + 3))); 
recC[nu_, mu_, k_, gsq_] = (-gsq)*(((mu + nu + 2*k + 1)*(mu + nu + 2*k + 2))/
     ((2*nu + 4*k + 3)*(2*nu + 4*k + 5))); 
cfA[nu_, mu_, k_, gsq_] := recC[nu, mu, k - 1, gsq]*recA[nu, mu, k, gsq]
cfB[nu_, mu_, k_, gsq_] := recB[nu, mu, k, gsq]

$matrixSize = 20; 

mat[n_, m_, gsq_, M_] := 
  Module[{kmin}, kmin = Round[Max[(m - n + Mod[n - m, 2])/2, -M]]; 
     Table[Switch[j - k, -1, recC[n, m, j, gsq], 0, recB[n, m, j, gsq], 1, 
       recA[n, m, j, gsq], _, 0], {j, kmin, M}, {k, kmin, M}]] /; 
   intQ[m, n] && 0 <= m <= n

mat[n_, m_, gsq_, M_] := Module[{kmin, kmax}, 
    kmin = Round[Max[(-m - n + Mod[n - m, 2])/2, -M]]; 
     kmax = Round[Min[(m - n - 2 + Mod[n - m, 2])/2, M]]; 
     Table[Switch[j - k, -1, recC[n, m, j, gsq], 0, recB[n, m, j, gsq], 1, 
       recA[n, m, j, gsq], _, 0], {j, kmin, kmax}, {k, kmin, kmax}]] /; 
   intQ[m, n] && m >= 0 && -m - 1 < n < m

mat[nu_, mu_, gsq_, M_] := Table[Switch[j - k, -1, recC[nu, mu, j, gsq], 0, 
    recB[nu, mu, j, gsq], 1, recA[nu, mu, j, gsq], _, 0], {j, -M, M}, 
   {k, -M, M}]

matsize[abs_] := Which[abs <= 100, 20, 100 < abs < 1000, 50, 1000 <= abs, 100]

startingEV[n_, m_, gsq_] := 
  Module[{M}, M = matsize[Abs[gsq]]; Sort[Eigenvalues[N[mat[n, m, gsq, M]]]][[
      -M - 1]]] /; intQ[m, n] && 0 <= m <= n

startingEV[n_, m_, gsq_] := 
  Module[{M, pos}, M = matsize[Abs[gsq]]; 
     pos = If[2*M + 1 < m, -M - 1, n + 1]; 
     Sort[Eigenvalues[N[mat[n, m, gsq, M]]]][[pos]]] /; 
   intQ[m, n] && 0 <= n < m

startingEV[n_, m_, gsq_] := startingEV[-n - 1, m, gsq] /; 
   intQ[m, n] && m >= 0 && n < 0
startingEV[n_, m_, gsq_] := startingEV[n, -m, gsq] /; intQ[m, n] && m < 0

startingEV[nu_, mu_, gsq_] := Module[{M, prec, nu1, mu1, gsq1, pos}, 
   M = matsize[Abs[gsq]]; prec = If[Abs[Mod[nu, 1] - 0.5] < 0.1, 17, 16]; 
    {nu1, mu1, gsq1} = SetPrecision[{nu, mu, gsq}, prec]; 
    pos = First[Flatten[Position[Sort[Table[(nu + 2*k)*(nu + 2*k + 1), 
          {k, -M, M}]], nu*(nu + 1)]]]; 
    Sort[Eigenvalues[N[mat[nu1, mu1, gsq1, M], prec]]][[pos]]]

startingEV[nu_, mu_, gsq_, ev_] := Module[{M, prec, nu1, mu1, gsq1}, 
   M = matsize[Abs[gsq]]; prec = If[Abs[Mod[nu, 1] - 0.5] < 0.1, 17, 16]; 
    {nu1, mu1, gsq1} = SetPrecision[{nu, mu, gsq}, prec]; 
    closestElt[Sort[Eigenvalues[N[mat[nu1, mu1, gsq1, M], prec]]], ev]]

U1deltaU1[nu_, mu_, gsq_, ev_] := 
   Module[{a, b, cf1, cf2, kmax, Nup, dU, res}, 
   a[k_] = cfA[nu, mu, k, gsq]; b[k_] = cfB[nu, mu, k, gsq] - ev; cf1 = 1; 
    cf2 = 0; kmax = 0; While[cf1 != cf2, kmax += 10; cf1 = cf2; 
      Nup[-kmax] = b[-kmax - 1]; Do[Nup[k + 1] = b[k] - a[k]/Nup[k], 
       {k, -kmax, 0}]; cf2 = Nup[1]]; dU = 1; Do[dU = 1 + dU*(a[k]/Nup[k]^2), 
     {k, -kmax, 0}]; res=Nup[1]; Clear[Nup]; {res, dU}]
U2deltaU2[nu_, mu_, gsq_, ev_] := 
   Module[{a, b, cf1, cf2, kmax, Ndown, dU, res}, 
   a[k_] = cfA[nu, mu, k, gsq]; b[k_] = cfB[nu, mu, k, gsq] - ev; cf1 = 1; 
    cf2 = 0; kmax = 0; While[cf1 != cf2, kmax += 10; cf1 = cf2; 
      Ndown[kmax] = a[kmax]/b[kmax]; 
      Do[Ndown[k] = a[k]/(b[k] - Ndown[k + 1]), {k, kmax - 1, 1, -1}]; 
      cf2 = Ndown[1]]; dU = 1; 
    Do[dU = If[a[k] == 0, 1, 1 + dU*(Ndown[k]^2/a[k])], {k, kmax, 1, -1}]; 
    dU -= 1; res=-Ndown[1]; Clear[Ndown]; {res, dU}]

bouwkamp[nu_, mu_, gsq_, ev_, prec_] := bouwkamp[nu, mu, gsq, ev, prec] =
  Module[{prec1, nu1, mu1, gsq1, ev1, ev2, U1, dU1, U2, dU2}, 
   prec1 = precRound[prec + $precIncInit]; While[True, {nu1, mu1, gsq1, ev1} = 
       SetPrecision[{nu, mu, gsq, ev}, prec1]; ev2 = 0; 
      While[ev1 != ev2, ev2 = ev1; {U1, dU1} = U1deltaU1[nu1, mu1, gsq1, 
          ev1]; {U2, dU2} = U2deltaU2[nu1, mu1, gsq1, ev1]; 
        ev1 += (U1 + U2)/(dU1 + dU2)]; If[Precision[ev1] >= prec, Break[]]; 
      prec1 += $precInc]; N[ev1, prec]]

numEV[nu_, mu_, gsq_, prec_] := bouwkamp[nu, mu, gsq, startingEV[nu, mu, gsq], prec]
numEV[nu_, mu_, gsq_, ev_, prec_] := bouwkamp[nu, mu, gsq, startingEV[nu, mu, gsq, ev], prec]

SpheroidalEigenvalue[nu_, mu_, g_] := N[nu*(nu + 1), Precision[{mu, g}]] /; 
   g == 0

SpheroidalEigenvalue[n_, m_, g_] := N[0, Precision[{m, n}]] /; 
   eltQ[n, {-1, 0}] && eltQ[m, {-1, 1}]
SpheroidalEigenvalue[n_, m_, g_] := 
  N[1 - Sqrt[1 - 4*g^2], Precision[{m, n}]] /; eltQ[n, {-1, 0}] && 
    eltQ[m, {-2, 2}]
SpheroidalEigenvalue[n_, m_, g_] := 
  N[1 + Sqrt[1 - 4*g^2], Precision[{m, n}]] /; eltQ[n, {-2, 1}] && 
    eltQ[m, {-2, 2}]
SpheroidalEigenvalue[n_, m_, g_] := 
  N[8/3 + ((1 + I*Sqrt[3])*(-28 + 48*g^2))/
      (12*(10 + 36*g^2 + 3*Sqrt[3]*Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^
        (1/3)) - (1 - I*Sqrt[3])*
      ((10 + 36*g^2 + 3*Sqrt[3]*Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^(1/3)/
       3), Precision[{m, n}]] /; eltQ[n, {-1, 0}] && eltQ[m, {-3, 3}]
SpheroidalEigenvalue[n_, m_, g_] := 
  N[8/3 + ((1 - I*Sqrt[3])*(-28 + 48*g^2))/
      (12*(10 + 36*g^2 + 3*Sqrt[3]*Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^
        (1/3)) - (1 + I*Sqrt[3])*
      ((10 + 36*g^2 + 3*Sqrt[3]*Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^(1/3)/
       3), Precision[{m, n}]] /; eltQ[n, {-2, 1}] && eltQ[m, {-3, 3}]
SpheroidalEigenvalue[n_, m_, g_] := 
  N[8/3 - (-28 + 48*g^2)/(6*(10 + 36*g^2 + 3*Sqrt[3]*
          Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^(1/3)) + 
     (2/3)*(10 + 36*g^2 + 3*Sqrt[3]*Sqrt[-9 + 92*g^2 - 64*g^4 + 64*g^6])^
       (1/3), Precision[{m, n}]] /; eltQ[n, {-3, 2}] && eltQ[m, {-3, 3}]

SpheroidalEigenvalue[n_, m_, g_] := N[0, Precision[{n, m, g}]] /; 
   natQ[n - 1] && eltQ[m, {-1, 1}] && g == n*(Pi/2)
SpheroidalEigenvalue[n_, m_, g_] := N[0, Precision[{n, m, g}]] /; 
   natQ[-n - 2] && eltQ[m, {-1, 1}] && g == (-n - 1)*(Pi/2)

SpheroidalEigenvalue[nu_, mu_, g_] := 
  N[MathieuCharacteristicA[nu + 1/2, g^2/4] - 1/4 - g^2/2,
    Precision[{nu, mu, g}]] /; eltQ[mu, {-2^(-1), 1/2}]

notHalfInt[nu_] := (Mod[nu, 1] != 1/2)


SpheroidalEigenvalue[nu_, mu_, g_] := 
  numEV[nu, mu, g^2, Precision[{nu, mu, g}]] /;
   numQ[nu, mu, g] && Precision[{nu, mu, g}] < Infinity && notHalfInt[nu] 
SpheroidalEigenvalue[nu_, mu_, g_, ev_] := 
  numEV[nu, mu, g^2, ev, Precision[{nu, mu, g, ev}]] /; 
   numQ[nu, mu, g, ev] && Precision[{nu, mu, g, ev}] < Infinity && 
    notHalfInt[nu]

unnormcoeff[nu_, mu_, gsq_, prec_, M_] := 
  Module[{prec1, ev, nu1, mu1, gsq1, A, B, C, up, down, list}, 
   prec1 = precRound[prec + $precIncInit];
   ev = N[SpheroidalEigenvalue[nu, mu, Sqrt[gsq]], prec1]; 
    While[True, {nu1, mu1, gsq1} = SetPrecision[{nu, mu, gsq}, prec1]; 
      A[k_] = recA[nu1, mu1, k, gsq1] /. zerorule; 
      B[k_] = recB[nu1, mu1, k, gsq1] - ev /. zerorule; 
      C[k_] = recC[nu1, mu1, k, gsq1] /. zerorule; up = down = {0}; 
      Do[down = (-A[k]/(B[k] + C[k]*First[down]))*Prepend[down, 1], 
       {k, M, 1, -1}]; Do[up = (-C[k]/(B[k] + A[k]*Last[up]))*Append[up, 1], 
       {k, -M, -1}]; list = Take[Join[up, {1}, down], {2, -2}] /. zerorule; 
      If[Precision[list] < prec || MemberQ[list, Indeterminate], 
       prec1 += $precInc; ev = bouwkamp[nu1, mu1, gsq1, ev, prec1], Break[]]]; 
    N[list, prec]]

pochterm[nu_, mu_, k_] := If[Pochhammer[nu - mu + 1, k] == 0 || 
    Pochhammer[nu + mu + 1, k] == ComplexInfinity, 0, 
   ((2*nu + 1)/(2*nu + 2*k + 1))*(Pochhammer[nu + mu + 1, k]/
     Pochhammer[nu - mu + 1, k])]

angcoeff[nu_, mu_, gsq_, prec_, M_] := angcoeff[nu, mu, gsq, prec, M] =
  (Module[{nu1, mu1, list, poch, sgn}, 
   {nu1, mu1} = SetPrecision[{nu, mu}, prec]; 
    list = unnormcoeff[nu, mu, gsq, prec, M]; 
    poch = Table[pochterm[nu1, mu1, 2*k], {k, -M, M}]; 
    sgn = If[list[[-M - 1]] == 0, 1, Sign[list[[-M - 1]]]];
    list*(sgn/(list^2 . poch)^(1/2))])

unnormcoeff2[nu_, mu_, gsq_, prec_, M_] := 
  (Module[{del, list, prec1, ev, kmin, kmax, nu1, mu1, gsq1, A, B, C, list2, 
     C2}, del = Round[Mod[mu + nu, 2]]; 
     prec1 = precRound[prec + $precIncInit];
     list = Take[unnormcoeff[nu, mu, gsq, prec1, M], 
       Max[-2*M - 1, -M - 1 + (-mu - nu + del)/2]];
     ev = SpheroidalEigenvalue[nu, mu, SetPrecision[Sqrt[gsq],prec1]];
     kmin = (Round[-mu - nu - 2] + del)/2; kmax = kmin - M; 
     While[True, {nu1, mu1, gsq1} = SetPrecision[{nu, mu, gsq}, prec1]; 
       A[k_] = recA[nu1, mu1, k, gsq1] /. zerorule; 
       B[k_] = recB[nu1, mu1, k, gsq1] - ev /. zerorule; 
       C[k_] = recC[nu1, mu1, k, gsq1] /. zerorule; list2 = {0}; 
       Do[list2 = (-C[k]/(B[k] + A[k]*First[list2]))*Prepend[list2, 1], 
        {k, kmax, kmin - 1}]; C2 = (-1)^del*(gsq1/((2*mu1 - 2*del - 1)*
           (2*mu1 - 2*del + 1))); list2 = 
        (-C2/(B[kmin] + A[kmin]*First[list2]))*Prepend[list2, 1]*First[list]; 
       list2 = N[Drop[list2, -1] /. zerorule, prec]; 
       If[Precision[list2] < prec || MemberQ[list2, Indeterminate], 
        prec1 += $precInc; ev = bouwkamp[nu, mu, gsq, ev, prec1];
      list = Take[unnormcoeff[nu, mu, gsq, prec1, M], 
                  Max[-2*M - 1, -M - 1 + (-mu - nu + del)/2]], Break[]]]; 
     {list, list2}]) /; natQ[mu + nu]

angcoeff2[nu_, mu_, gsq_, prec_, M_] := 
  (angcoeff2[nu, mu, gsq, prec, M] = 
    Module[{nu1, mu1, list, list2, poch, sgn}, 
    {nu1, mu1} = SetPrecision[{nu, mu}, prec + $precIncInitSmall]; 
     {list, list2} = unnormcoeff2[nu, mu, gsq, prec + $precIncInitSmall, M]; 
     poch = Table[pochterm[nu1, mu1, 2*k], 
       {k, Max[-M, Round[-mu1 - nu1 + Mod[mu1 + nu1, 2]]/2], M}]; 
     sgn = If[list[[-M - 1]] == 0, 1, Sign[list[[-M - 1]]]]; 
     {list, list2}*(sgn/(list^2 . poch)^(1/2))]) /; natQ[mu + nu]

normfactor[nu_, mu_, gsq_, prec_] := normfactor[nu, mu, gsq, prec] =
Module[{prec1, M, list, A, tail}, prec1 = precRound[prec + $precIncInitSmall]; M = $seriesLength; 
     While[True, list = Table[(-1)^k, {k, -M, M}]*angcoeff[nu, mu, gsq, 
          prec1, M]; A = Plus @@ list; 
       tail = Max[Abs[{First[list], Last[list]}]]; If[Precision[A] < prec, 
        prec1 += $precInc]; If[tail/Abs[A] > 10^(-prec), M += $seriesInc]; 
       If[Precision[A] >= prec && tail/Abs[A] < 10^(-prec), Break[]]]; A]

SpheroidalRadialFactor[nu_, mu_, g_] := 
  normfactor[nu, mu, g^2, precRound[Precision[{nu, mu, g}]]] /; 
   numQ[nu, mu, g] && Precision[{nu, mu, g}] < Infinity

$seriesLength = 40;

numericP[n_, m_, type_, g_, prec_, z_, legTable_] := 
  (Module[{M, prec1, n1, m1, kmin, kmax, 
      list, testlist, testvalue, tail, sgnlist, ltable},
      M = $seriesLength;
      prec1 = prec;
      {n1, m1} = Round[{n, m}]; 
      While[True, kmin = Max[-M, (-n1 + Abs[m1] + Mod[m1 + n1, 2])/2]; 
        kmax = M; sgnlist = Table[(-1)^k, {k, kmin, kmax}];
        list = Take[angcoeff[n1, m1, g^2, precRound[prec1 + $precIncInitSmall], M], kmin - kmax - 1];
        ltable = legTable[n1, m1, type, kmin, kmax, SetPrecision[z, prec1]];
        testlist = list*sgnlist*ltable;
        testvalue = Plus @@ testlist; 
        tail = Max[Abs[{If[-M > (-n1 + Abs[m1] + Mod[m1 + n1, 2])/2, 
             First[testlist], 0], Last[testlist]}]]; 
        If[Precision[testvalue] < prec, prec1 += $precInc]; 
        If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
        If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
         Break[]]]; N[testvalue, prec]])

numericPS[n_, m_, type_, g_, prec_, z_] :=
  numericP[n, m, type, g, prec, z, legPTableNoCheck] /; intQ[m, n] && n >= Abs[m]


numericPS[n_, m_, type_, g_, prec_, z_] := numericPS[-n - 1, m, type, g, prec, z] /; 
   intQ[m, n] && -n - 1 >= Abs[m]

numericPSPrime[n_, m_, type_, g_, prec_, z_] :=
  numericP[n, m, type, g, prec, z, legPPrimeTable] /; intQ[m, n] && n >= Abs[m]

numericPSPrime[n_, m_, type_, g_, prec_, z_] := 
  numericPSPrime[-n - 1, m, type, g, prec, z] /; intQ[m, n] && -n - 1 >= Abs[m]

numericQ[n_, m_, type_, g_, prec_, z_, legTable1_, legTable2_] :=
   (Module[{M, prec1, del, list, 
      list2, kmin1, kmax1, kmin2, kmax2, ztest, testlist, testlist2, 
      testvalue, tail, sgnlist, sgnlist2, ltable1, ltable2},
     M = $seriesLength; prec1 = precRound[prec + $precIncInitSmall];
     del = Round[Mod[m + n, 2]];
      While[True,
       {list, list2} = angcoeff2[n, m, g^2, prec1, M];
        kmin1 = Max[-M, (-m - n + del)/2]; kmax1 = M;
        kmin2 = 0; kmax2 = Length[list2] - 1;
        sgnlist = Table[(-1)^k, {k, kmin1, kmax1}];
        sgnlist2 = (-1)^((m + n - del + 2)/2)*Table[(-1)^k, {k, kmin2, kmax2}];
        ztest = SetPrecision[z, prec1];
        ltable1 = legTable1[n, m, type, kmin1, kmax1, ztest]; 
        ltable2 = legTable2[m - del + 1, m, type, kmin2, kmax2, ztest]; 
        testlist = list*sgnlist*ltable1;
        testlist2 = list2*sgnlist2*ltable2;
        testvalue = Plus @@ Join[testlist, testlist2];
        tail = Max[Abs[{Last[testlist], Last[testlist2]}]]; 
        If[Precision[testvalue] < prec, prec1 += $precInc]; 
        If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
        If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec,
         Break[]]]; N[testvalue, prec]]) /; intQ[m, n] && n >= -m

numericQS[n_, m_, type_, g_, prec_, z_] :=
   numericQ[n, m, type, g, prec, z, legQTable, legPTable] /; intQ[m, n] && n >= -m

numericQSPrime[n_, m_, type_, g_, prec_, z_] :=
   numericQ[n, m, type, g, prec, z, legQPrimeTable, legPPrimeTable] /; intQ[m, n] && n >= -m

numericQ[nu_, mu_, type_, g_, prec_, z_, legTable1_, legTable2_, legTable3_] :=
   (Module[{M, prec1, del, nu1, mu1, list, 
      list2, kmin1, kmax1, kmin2, kmax2, ztest, testlist, testlist2, 
      testvalue, tail, sgnlist, sgnlist2, legtable1, legtable2, legtable3},
     M = $seriesLength; prec1 = precRound[prec + $precIncInitSmall];
     del = Round[Mod[mu + nu, 2]]; 
      While[True, {nu1, mu1} = SetPrecision[{nu, mu}, prec1]; 
        {list, list2} = angcoeff2[nu, mu, g^2, prec1, M]; 
        kmin1 = Max[-M, (-Round[mu + nu] + del)/2]; kmax1 = M; kmin2 = 0;
        kmax2 = Length[list2] - 1;
        sgnlist = Table[(-1)^k, {k, kmin1, kmax1}];
        sgnlist2 = Table[(-1)^k, {k, kmin2, kmax2}];
        ztest = SetPrecision[z, prec1];
        ltable1 = legTable1[nu1, mu1, type, kmin1, kmax1, ztest];
        ltable2 = legTable2[mu1 - del + 1, mu1, type, kmin2, kmax2, ztest];
        ltable3 = legTable3[mu1 - del + 1, mu1, type, kmin2, kmax2, ztest];
        testlist = list*sgnlist*ltable1;
        testlist2 = (-1)^del*list2*sgnlist2*
          (Exp[I*mu1*Pi]*Cos[nu1*Pi]*ltable2 + (Sin[Pi*(nu1 - mu1)]/Pi)*ltable3);
        testvalue = Plus @@ Join[testlist, testlist2]; 
        tail = Max[Abs[{Last[testlist], Last[testlist2]}]]; 
        If[Precision[testvalue] < prec, prec1 += $precInc]; 
        If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
        If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
         Break[]]]; N[testvalue, prec]]) /; natQ[mu + nu] && !intQ[mu, nu]

numericQS[nu_, mu_, type_, g_, prec_, z_] :=
    numericQ[nu, mu, type, g, prec, z, legQTable, legPTable, legQTable]  /; natQ[mu + nu] && !intQ[mu, nu]

numericQSPrime[nu_, mu_, type_, g_, prec_, z_] :=
    numericQ[nu, mu, type, g, prec, z, legQPrimeTable, legPPrimeTable, legQTable]  /; natQ[mu + nu] && !intQ[mu, nu]

numericPQ[nu_, mu_, type_, g_, prec_, z_, legTable_] := (
    (Module[{M, prec1, list, testlist, testvalue, tail, ltable}, 
      M = $seriesLength; prec1 = precRound[prec + $precIncInitSmall]; While[True,
      list = Table[(-1)^k, {k, -M, M}]*angcoeff[nu, mu, g^2, prec1, M];
      ltable = legTable[nu, mu, type, -M, M, SetPrecision[z, prec1]];
      testlist = list*ltable;
      testvalue = Plus @@ testlist; 
       tail = Max[Abs[{First[testlist], Last[testlist]}]]; 
       If[Precision[testvalue] < prec, prec1 += $precInc]; 
       If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
       If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
        Break[]]]; N[testvalue, prec]]))

numericPS[nu_, mu_, type_, g_, prec_, z_] := numericPQ[nu, mu, type, g, prec, z, legPTable]
numericQS[nu_, mu_, type_, g_, prec_, z_] := numericPQ[nu, mu, type, g, prec, z, legQTable]
numericPSPrime[nu_, mu_, type_, g_, prec_, z_] := numericPQ[nu, mu, type, g, prec, z, legPPrimeTable]
numericQSPrime[nu_, mu_, type_, g_, prec_, z_] := numericPQ[nu, mu, type, g, prec, z, legQPrimeTable]

numericS[n_, m_, g_, prec_, z_, besTable_] := (
      Module[{M, prec1, kmin, kmax, list, testlist, testvalue, tail, A, btable},
      M = $seriesLength; prec1 = prec;
      While[True, kmin = Max[-M, (-Abs[m] - n + Mod[n + m, 2])/2]; kmax = M; 
        list = Take[angcoeff[n, -m, g^2, precRound[prec1 + $precIncInitSmall], M], kmin - kmax - 1];
        A = Table[(-1)^k, {k, kmin, kmax}] . list;
        btable = besTable[n, kmin, kmax, SetPrecision[g*z, prec1]]; 
        testlist = (1/A)*list*btable;
        testvalue = Plus @@ testlist; 
        tail = Max[Abs[{If[-M > (-Abs[m] - n + Mod[n + m, 2])/2, 
             First[testlist], 0], Last[testlist]}]]; 
        If[Precision[testvalue] < prec, prec1 += $precInc]; 
        If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
        If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
         Break[]]]; If[m != 0, N[(1 - 1/z^2)^(m/2)*testvalue, prec], 
        N[testvalue, prec]]]) /; intQ[m, n] && n >= Abs[m]

numericSPrime[n_, m_, g_, prec_, z_, besTable_, besPrimeTable_] :=
     (Module[{M, prec1, kmin, kmax, list, testlist, testvalue, tail, A, btable},
      M = $seriesLength; 
      prec1 = precRound[prec + $precIncInitSmall];
      While[True, kmin = Max[-M, (-Abs[m] - n + Mod[n + m, 2])/2]; kmax = M;
        list = Take[angcoeff[n, -m, g^2, prec1, M], kmin - kmax - 1];
        btable = besPrimeTable[n, kmin, kmax, SetPrecision[g*z, prec1]];
        A = Table[(-1)^k, {k, kmin, kmax}] . list; 
        testlist = (1/A)*list*btable;
        testvalue = Plus @@ testlist; 
        tail = Max[Abs[{If[-M > (-Abs[m] - n + Mod[n + m, 2])/2, 
             First[testlist], 0], Last[testlist]}]]; 
        If[Precision[testvalue] < prec, prec1 += $precInc]; 
        If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
        If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
         Break[]]]; If[m != 0,
        N[(g*(1 - 1/z^2)^(m/2)*testvalue + m*(1 - 1/z^2)^(m/2 - 1)*
            (1/A)*(list . besTable[n, kmin, kmax, SetPrecision[g*z, prec1]]/z^3)), prec], 
        N[g*testvalue, prec]]]) /; intQ[m, n] && n >= Abs[m]

numericS[nu_, mu_, g_, prec_, z_, besTable_] := (
    Module[{M, prec1, list, testlist, testvalue, tail, A, btable}, 
    M = $seriesLength; prec1 = precRound[prec + $precIncInitSmall]; While[True, 
      A = normfactor[nu, -mu, g^2, prec1];
       list = angcoeff[nu, -mu, g^2, prec1, M]/A;
       btable = besTable[nu, -M, M, SetPrecision[g*z, prec1]];
       testlist = list*btable;
       testvalue = Plus @@ testlist; 
       tail = Max[Abs[{First[testlist], Last[testlist]}]]; 
       If[Precision[testvalue] < prec, prec1 += $precInc]; 
       If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
       If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
        Break[]]]; If[mu != 0,
       N[(1 - 1/z^2)^(mu/2)*testvalue, prec],
       N[testvalue, prec]]])

numericSPrime[nu_, mu_, g_, prec_, z_, besTable_, besPrimeTable_] := 
    Module[{M, prec1, list, testlist, testvalue, tail, A, btable}, 
    M = $seriesLength; prec1 = precRound[prec + $precIncInitSmall]; While[True, 
      A = normfactor[nu, -mu, g^2, prec1];
       list = angcoeff[nu, -mu, g^2, prec1, M]/A;
       btable = besPrimeTable[nu, -M, M, SetPrecision[g*z, prec1]];
       testlist = list*btable;
       testvalue = Plus @@ testlist; 
       tail = Max[Abs[{First[testlist], Last[testlist]}]]; 
       If[Precision[testvalue] < prec, prec1 += $precInc]; 
       If[tail > Abs[testvalue]/10^prec, M += $seriesInc]; 
       If[Precision[testvalue] >= prec && tail <= Abs[testvalue]/10^prec, 
        Break[]]]; If[mu != 0, 
       N[g*(1 - 1/z^2)^(mu/2)*testvalue + mu*(1 - 1/z^2)^(mu/2 - 1)*
          (list . besTable[nu, -M, M, SetPrecision[g*z, prec1]]/z^3), 
        prec], N[g*testvalue, prec]]]

numericS1[nu_, mu_, g_, prec_, z_] := numericS[nu, mu, g, prec, z, besJTable]
numericS2[nu_, mu_, g_, prec_, z_] := numericS[nu, mu, g, prec, z, besYTable]
numericS3[nu_, mu_, g_, prec_, z_] := numericS[nu, mu, g, prec, z, besKTable]

numericS1Prime[nu_, mu_, g_, prec_, z_] := numericSPrime[nu, mu, g, prec, z, besJTable, besJPrimeTable]
numericS2Prime[nu_, mu_, g_, prec_, z_] := numericSPrime[nu, mu, g, prec, z, besYTable, besYPrimeTable]
numericS3Prime[nu_, mu_, g_, prec_, z_] := numericSPrime[nu, mu, g, prec, z, besKTable, besKPrimeTable]

SpheroidalPS[nu_, mu_, type_:1, z_, g_] := LegendreP[nu, mu, type, z] /; 
   g == 0 && MemberQ[{1, 2}, type]
SpheroidalQS[nu_, mu_, type_:1, z_, g_] := LegendreQ[nu, mu, type+1, z] /; 
   g == 0 && MemberQ[{1, 2}, type]
SpheroidalPSPrime[nu_, mu_, type_:1, z_, g_] := 
  Module[{zz}, D[LegendreP[nu, mu, type, zz], zz] /. zz -> z] /; 
   g == 0 && MemberQ[{1, 2}, type]
SpheroidalQSPrime[nu_, mu_, type_:1, z_, g_] := 
  Module[{zz}, D[LegendreQ[nu, mu, type+1, zz], zz] /. zz -> z] /; 
   g == 0 && MemberQ[{1, 2}, type]

angNormOmega[n_] := 
  Sqrt[(n*(n + 1))/((2*n + 1)*(EulerGamma - CosIntegral[2*n*Pi] + 
      Log[2*n*Pi]))]
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  N[(-1)^((n + 1)/2)*angNormOmega[n]*(Cos[(n*(Pi/2))*z]/Sqrt[1 - z^2]), 
    Precision[{m, g}]] /; natQ[n - 1] && OddQ[Round[n]] && m == 1 && 
    g == n*(Pi/2) &&  !eltQ[z, {-1, 1}] && type === 1
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  N[(-1)^(n/2)*angNormOmega[n]*(Sin[(n*(Pi/2))*z]/Sqrt[1 - z^2]), 
    Precision[{m, g}]] /; natQ[n - 1] && EvenQ[Round[n]] && m == 1 && 
    g == n*(Pi/2) &&  !eltQ[z, {-1, 1}] && type === 1
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  0 /; natQ[n - 1] && m == 1 && g == n*(Pi/2) && eltQ[z, {-1, 1}] && 
    type === 1
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  (-(n - 1)!/(n + 1)!)*SpheroidalPS[n, -m, type, z, g] /; 
   natQ[n - 1] && m == -1 && g == n*(Pi/2) && type === 1
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  SpheroidalPS[-n - 1, m, type, z, g] /; natQ[-n - 2] && eltQ[m, {-1, 1}] && 
    g == (-n - 1)*(Pi/2) && type === 1
SpheroidalPS[n_, m_, 2, z_, g_] := 
  ((1 - z)^(m/2)/(z - 1)^(m/2))*SpheroidalPS[n, m, 1, z, g] /; 
   (natQ[n - 1] && g == n*(Pi/2) || natQ[-n - 2] && g == (-n - 1)*(Pi/2)) && 
    eltQ[m, {-1, 1}]
SpheroidalPSPrime[n_, m_, type_:1, z_, g_] := 
  Module[{zz}, D[SpheroidalPS[n, m, type, zz, g], zz] /. zz -> z] /; 
   (natQ[n - 1] && g == n*(Pi/2) || natQ[-n - 2] && g == (-n - 1)*(Pi/2)) && 
    eltQ[m, {-1, 1}] && MemberQ[{1, 2}, type]
SpheroidalQS[n_, m_, type_:1, z_, g_] := 
  N[(-1)^((n + 1)/2)*angNormOmega[n]*(Sin[(n*(Pi/2))*z]/Sqrt[1 - z^2]), 
    Precision[{m, g}]] /; natQ[n - 1] && OddQ[Round[n]] && m == 1 && 
    g == n*(Pi/2) &&  !eltQ[z, {-1, 1}] && type === 1
SpheroidalQS[n_, m_, type_:1, z_, g_] := 
  N[(-1)^((n + 2)/2)*angNormOmega[n]*(Cos[(n*(Pi/2))*z]/Sqrt[1 - z^2]), 
    Precision[{m, g}]] /; natQ[n - 1] && EvenQ[Round[n]] && m == 1 && 
    g == n*(Pi/2) &&  !eltQ[z, {-1, 1}] && type === 1
SpheroidalQS[n_, m_, type_:1, z_, g_] := ComplexInfinity /; 
   natQ[n - 1] && m == 1 && g == n*(Pi/2) && eltQ[z, {-1, 1}] && type === 1
SpheroidalQS[n_, m_, type_:1, z_, g_] := 
  (-(n - 1)!/(n + 1)!)*SpheroidalQS[n, -m, type, z, g] /; 
   natQ[n - 1] && m == -1 && g == n*(Pi/2) && type === 1
SpheroidalQS[n_, m_, 2, z_, g_] := (-1)^m*((z - 1)^(m/2)/(1 - z)^(m/2))*
    (SpheroidalQS[n, m, 1, z, g] + (Pi/2)*(Sqrt[1 - z]/Sqrt[z - 1])*
      SpheroidalPS[n, m, 1, z, g]) /; 
   (natQ[n - 1] && g == n*(Pi/2) || natQ[-n - 2] && g == (-n - 1)*(Pi/2)) && 
    eltQ[m, {-1, 1}]
SpheroidalQSPrime[n_, m_, type_:1, z_, g_] := 
  Module[{zz}, D[SpheroidalQS[n, m, type, zz, g], zz] /. zz -> z] /; 
   (natQ[n - 1] && g == n*(Pi/2) || natQ[-n - 2] && g == (-n - 1)*(Pi/2)) && 
    eltQ[m, {-1, 1}]

SpheroidalS1[n_, m_, z_, g_] := (1/(n*(Pi/2)*z*Sqrt[1 - 1/z^2]))*
    Sin[n*(Pi/2)*(z - 1)] /; eltQ[m, {-1, 1}] && natQ[n - 1] && g == n*(Pi/2)
SpheroidalS2[n_, m_, z_, g_] := (-(n*(Pi/2)*z*Sqrt[1 - 1/z^2])^(-1))*
    Cos[n*(Pi/2)*(z - 1)] /; eltQ[m, {-1, 1}] && natQ[n - 1] && g == n*(Pi/2)
SpheroidalS1[n_, m_, z_, g_] := (-1)^n*SpheroidalS2[-n - 1, m, z, g] /; 
   eltQ[m, {-1, 1}] && natQ[-n - 2] && g == (-n - 1)*(Pi/2)
SpheroidalS2[n_, m_, z_, g_] := (-1)^(n + 1)*SpheroidalS1[-n - 1, m, z, g] /; 
   eltQ[m, {-1, 1}] && natQ[-n - 2] && g == (-n - 1)*(Pi/2)
SpheroidalS1Prime[n_, m_, z_, g_] := 
  Module[{zz}, D[SpheroidalS1[n, m, zz, g], zz] /. zz -> z] /; 
   eltQ[m, {-1, 1}] && (natQ[n - 1] && g == n*(Pi/2) || 
     natQ[-n - 2] && g == (-n - 1)*(Pi/2))
SpheroidalS2Prime[n_, m_, z_, g_] := 
  Module[{zz}, D[SpheroidalS2[n, m, zz, g], zz] /. zz -> z] /; 
   eltQ[m, {-1, 1}] && (natQ[n - 1] && g == n*(Pi/2) || 
     natQ[-n - 2] && g == (-n - 1)*(Pi/2))

SpheroidalPS[nu_, mu_, type_:1, z_, g_] := 
  N[Sqrt[2/Pi]*(1/(1 - z^2)^(1/4))*MathieuC[MathieuCharacteristicA[nu + 1/2, 
       g^2/4], g^2/4, ArcCos[z]], Precision[{mu, nu}]] /; 
   mu == 1/2 && notHalfInt[nu] && type === 1
SpheroidalQS[nu_, mu_, type_:1, z_, g_] := 
  N[(-Sqrt[Pi/2])*(1/(1 - z^2)^(1/4))*MathieuS[MathieuCharacteristicB[
       nu + 1/2, g^2/4], g^2/4, ArcCos[z]], Precision[{mu, nu}]] /; 
   mu == 1/2 && notHalfInt[nu] && type === 1
SpheroidalPS[nu_, mu_, type_:1, z_, g_] := 
  N[(1/Sqrt[Pi])*(1/(1 - z^2)^(1/4))*MathieuC[MathieuCharacteristicA[
       nu + 1/2, g^2/4], g^2/4, ArcCos[z]], Precision[{mu, nu}]] /; 
   mu == 1/2 && Mod[nu, 1] == 1/2 && type === 1
SpheroidalQS[nu_, mu_, type_:1, z_, g_] := 
  0 /; mu == 1/2 && Mod[nu, 1] == 1/2 && type === 1
SpheroidalPS[nu_, mu_, type_:1, z_, g_] := 
  (-Gamma[nu + 1/2]/Gamma[nu + 3/2])*(2/Pi)*SpheroidalQS[nu, -mu, type, z, 
     g] /; mu == -2^(-1)
SpheroidalQS[nu_, mu_, type_:1, z_, g_] := 
  (Gamma[nu + 1/2]/Gamma[nu + 3/2])*(Pi/2)*SpheroidalPS[nu, -mu, type, z, 
     g] /; mu == -2^(-1)
SpheroidalPS[nu_, mu_, 2, z_, g_] := 
  ((1 - z)^(mu/2)/(z - 1)^(mu/2))*SpheroidalPS[nu, mu, 1, z, g] /; mu^2 == 1/4
SpheroidalQS[nu_, mu_, 2, z_, g_] := 
  Exp[mu*Pi*I]*((z - 1)^(mu/2)/(1 - z)^(mu/2))*
    (SpheroidalQS[nu, mu, 1, z, g] + Pi*(Csc[mu*Pi]/2)*
      ((1 - z)^mu/(z - 1)^mu - Cos[mu*Pi])*SpheroidalPS[nu, mu, 1, z, g]) /; 
   mu^2 == 1/4
SpheroidalPSPrime[nu_, mu_, type_:1, z_, g_] := 
  Module[{zz}, D[SpheroidalPS[nu, mu, type, zz, g], zz] /. zz -> z] /; 
   mu^2 == 1/4
SpheroidalQSPrime[nu_, mu_, type_:1, z_, g_] := 
  Module[{zz}, D[SpheroidalQS[nu, mu, type, zz, g], zz] /. zz -> z] /; 
   mu^2 == 1/4

SpheroidalQS[nu_, mu_, type_:1, z_, g_] := ComplexInfinity /; 
   natQ[-mu - nu - 1]
SpheroidalQSPrime[nu_, mu_, type_:1, z_, g_] := 
  ComplexInfinity /; natQ[-mu - nu - 1]
SpheroidalPS[n_, m_, type_:1, z_, g_] := 
  0 /; intQ[m, n] && -Abs[m] - 1 < n < Abs[m]
SpheroidalPSPrime[n_, m_, type_:1, z_, g_] :=
  0 /; intQ[m, n] && -Abs[m] - 1 < n < Abs[m]


outsideRadius[z_, g_] := (Abs[z] > 2)

SpheroidalPS[nu_, mu_, type_:1, z_, g_] := 
  numericPS[nu, mu, type, g, precRound[Precision[{nu, mu, z, g}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, z, g}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] && MemberQ[{1, 2}, type]
SpheroidalPSPrime[nu_, mu_, type_:1, z_, g_] := 
  numericPSPrime[nu, mu, type, g, precRound[Precision[{nu, mu, z, g}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, z, g}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] && MemberQ[{1, 2}, type]
SpheroidalQS[nu_, mu_, type_:1, z_, g_] := 
  numericQS[nu, mu, type, g, precRound[Precision[{nu, mu, z, g}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, z, g}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] && MemberQ[{1, 2}, type]
SpheroidalQSPrime[nu_, mu_, type_:1, z_, g_] := 
  numericQSPrime[nu, mu, type, g, precRound[Precision[{nu, mu, z, g}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, z, g}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] && MemberQ[{1, 2}, type]

SpheroidalS1[nu_, mu_, z_, g_] := 
  numericS1[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] && g != 0
SpheroidalS2[nu_, mu_, z_, g_] := 
  numericS2[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu]
SpheroidalS3[nu_, mu_, z_, g_] := 
  numericS3[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu]
SpheroidalS1Prime[nu_, mu_, z_, g_] := 
  numericS1Prime[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] && g != 0
SpheroidalS2Prime[nu_, mu_, z_, g_] := 
  numericS2Prime[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu]
SpheroidalS3Prime[nu_, mu_, z_, g_] := 
  numericS3Prime[nu, mu, g, precRound[Precision[{nu, mu, g, z}]], z] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu]

SpheroidalPS[nu_, mu_, type_:1, z_, g_] := 
  ((z - 1)^(mu/2)/(1 - z)^(mu/2))*SpheroidalPS[nu, mu, 2, z, g] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] && type === 1
SpheroidalPSPrime[nu_, mu_, type_:1, z_, g_] := 
  ((z - 1)^(mu/2)/(1 - z)^(mu/2))*SpheroidalPSPrime[nu, mu, 2, z, g] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] && type === 1

SpheroidalQS[nu_, mu_, type_:1, z_, g_] := 
  ((1 - z)^(mu/2)/(z - 1)^(mu/2))*
    (Exp[(-I)*mu*Pi]*SpheroidalQS[nu, mu, 2, z, g] + 
     (Pi/2)*Csc[mu*Pi]*(Cos[mu*Pi]*((z - 1)^mu/(1 - z)^mu) - 1)*
      SpheroidalPS[nu, mu, 2, z, g]) /; numQ[nu, mu, g, z] && 
    Precision[{nu, mu, g, z}] < Infinity && outsideRadius[z,g] && 
    notHalfInt[nu] &&  !intQ[mu] && type === 1
SpheroidalQSPrime[nu_, mu_, type_:1, z_, g_] := 
  ((1 - z)^(mu/2)/(z - 1)^(mu/2))*
    (Exp[(-I)*mu*Pi]*SpheroidalQSPrime[nu, mu, 2, z, g] + 
     (Pi/2)*Csc[mu*Pi]*(Cos[mu*Pi]*((z - 1)^mu/(1 - z)^mu) - 1)*
      SpheroidalPSPrime[nu, mu, 2, z, g]) /; numQ[nu, mu, g, z] && 
    Precision[{nu, mu, g, z}] < Infinity && outsideRadius[z,g] && 
    notHalfInt[nu] &&  !intQ[mu] && type === 1

SpheroidalQS[nu_, m_, type_:1, z_, g_] := 
  (-1)^m*((1 - z)^(m/2)/(z - 1)^(m/2))*(SpheroidalQS[nu, m, 2, z, g] - 
     (Pi/2)*(Sqrt[1 - z]/Sqrt[z - 1])*SpheroidalPS[nu, m, 2, z, g]) /; 
   numQ[nu, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] && 
    notHalfInt[nu] && intQ[m] && type === 1
SpheroidalQSPrime[nu_, m_, type_:1, z_, g_] := 
  (-1)^m*((1 - z)^(m/2)/(z - 1)^(m/2))*(SpheroidalQSPrime[nu, m, 2, z, g] - 
     (Pi/2)*(Sqrt[1 - z]/Sqrt[z - 1])*SpheroidalPSPrime[nu, m, 2, z, g]) /; 
   numQ[nu, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] && 
    notHalfInt[nu] && intQ[m] && type === 1

SpheroidalQS[nu_, mu_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{nu, mu, g, z}]];
     (1/numericJF[-nu - 1, mu, g, prec])*Pi*Csc[(mu + nu)*Pi]*
      Exp[I*(mu - nu)*Pi]*(z - 1)^(mu/2)*(z + 1)^(mu/2)*
      ((g*z)^nu/((1 - 1/z^2)^(mu/2)*g^nu*z^(mu + nu)))*
      numericS1[-nu - 1, mu, g, prec, z]] /; numQ[nu, mu, g, z] && 
    Precision[{nu, mu, g, z}] < Infinity && outsideRadius[z,g] && 
    notHalfInt[nu] &&  !intQ[mu, nu]
SpheroidalQSPrime[nu_, mu_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{nu, mu, g, z}]]; 
     (1/numericJF[-nu - 1, mu, g, prec])*Pi*Csc[(mu + nu)*Pi]*
      Exp[I*(mu - nu)*Pi]*(z - 1)^(mu/2)*(z + 1)^(mu/2)*
      ((g*z)^nu/((1 - 1/z^2)^(mu/2)*g^nu*z^(mu + nu)))*
      numericS1Prime[-nu - 1, mu, g, prec, z]] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu]

SpheroidalPS[nu_, mu_, 2, z_, g_] := 
  (Sec[nu*Pi]/Pi)*Exp[(-mu)*Pi*I]*
    (Sin[(mu + nu)*Pi]*SpheroidalQS[nu, mu, 2, z, g] + 
     Sin[(mu - nu)*Pi]*SpheroidalQS[-nu - 1, mu, 2, z, g]) /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu]
SpheroidalPSPrime[nu_, mu_, 2, z_, g_] := 
  (Sec[nu*Pi]/Pi)*Exp[(-mu)*Pi*I]*
    (Sin[(mu + nu)*Pi]*SpheroidalQSPrime[nu, mu, 2, z, g] + 
     Sin[(mu - nu)*Pi]*SpheroidalQSPrime[-nu - 1, mu, 2, z, g]) /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity &&
    outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu]

radialPrefactor[m_,z_] := 1 /; EvenQ[m] || z==1
radialPrefactor[m_,z_] := (1 - 1/z^2)^(m/2)*(z^m/((z - 1)^(m/2)*(z + 1)^(m/2)))
angularPrefactor[m_,z_]:= 1/radialPrefactor[m,z]

SpheroidalQS[n_, m_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, g, z}]]; 
     (-1)^(m + 1)*g*numericJF[n, -m, g, prec]*normfactor[n, -m, g^2, prec]*
      normfactor[n, m, g^2, prec]*angularPrefactor[m,z]*numericS2[n, m, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] && 
    intQ[m, n]
SpheroidalQSPrime[n_, m_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, g, z}]]; 
     (-1)^(m + 1)*g*numericJF[n, -m, g, prec]*normfactor[n, -m, g^2, prec]*
      normfactor[n, m, g^2, prec]*angularPrefactor[m,z]*numericS2Prime[n, m, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] &&
    intQ[m, n]

SpheroidalPS[n_, m_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, g, z}]]; 
     (1/numericJF[n, m, g, prec])*angularPrefactor[m,z]*numericS1[n, m, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] && 
    intQ[m, n]
SpheroidalPSPrime[n_, m_, 2, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, g, z}]]; 
     (1/numericJF[n, m, g, prec])*angularPrefactor[m,z]*numericS1Prime[n, m, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && outsideRadius[z,g] && 
    intQ[m, n]


SpheroidalS1[n_, m_, z_, g_] := (-1)^(n) SpheroidalS2[-n-1, m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && n < 0

SpheroidalS1Prime[n_, m_, z_, g_] := (-1)^(n) SpheroidalS2Prime[-n-1, m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && n < 0

SpheroidalS2[n_, m_, z_, g_] := (-1)^(n+1) SpheroidalS1[-n-1, m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && n < 0

SpheroidalS2Prime[n_, m_, z_, g_] := (-1)^(n+1) SpheroidalS1Prime[-n-1, m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && n < 0

SpheroidalS1[nu_, mu_, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{nu, mu, z, g}]]; 
     numericJF[nu, mu, g, prec]*(Sin[(mu - nu)*Pi]/Pi)*Exp[(-I)*(mu + nu)*Pi]*
      (((1 - z^(-2))^(mu/2)*(g*z)^nu)/(g^nu*z^(nu - mu)*(z - 1)^(mu/2)*
        (z + 1)^(mu/2)))*numericQS[-nu - 1, mu, 2, g, prec, z]] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu]  && g != 0

SpheroidalS1Prime[nu_, mu_, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{nu, mu, z, g}]]; 
     numericJF[nu, mu, g, prec]*(Sin[(mu - nu)*Pi]/Pi)*Exp[(-I)*(mu + nu)*Pi]*
      (((1 - z^(-2))^(mu/2)*(g*z)^nu)/(g^nu*z^(nu - mu)*(z - 1)^(mu/2)*
        (z + 1)^(mu/2)))*numericQSPrime[-nu - 1, mu, 2, g, prec, z]] /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu] && g != 0

SpheroidalS2[nu_, mu_, z_, g_] := 
  (-Sec[nu*Pi])*(SpheroidalS1[-nu - 1, mu, z, g] + 
     Sin[nu*Pi]*SpheroidalS1[nu, mu, z, g]) /; numQ[nu, mu, g, z] && 
    Precision[{nu, mu, g, z}] < Infinity && !outsideRadius[z,g] && 
    notHalfInt[nu] &&  !intQ[mu, nu]

SpheroidalS2Prime[nu_, mu_, z_, g_] := 
  (-Sec[nu*Pi])*(SpheroidalS1Prime[-nu - 1, mu, z, g] + 
     Sin[nu*Pi]*SpheroidalS1Prime[nu, mu, z, g]) /; 
   numQ[nu, mu, g, z] && Precision[{nu, mu, g, z}] < Infinity && 
    !outsideRadius[z,g] && notHalfInt[nu] &&  !intQ[mu, nu]

SpheroidalS1[n_, m_, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, z, g}]]; 
     numericJF[n, m, g, prec]*radialPrefactor[m,z]*numericPS[n, m, 2, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && !outsideRadius[z,g] && 
    intQ[m, n]  && g != 0
SpheroidalS1Prime[n_, m_, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, z, g}]]; 
     numericJF[n, m, g, prec]*radialPrefactor[m, z]*numericPSPrime[n, m, 2, g,
      prec, z]]  /; numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && 
    !outsideRadius[z,g] && intQ[m, n] && g != 0

SpheroidalS2[n_, m_, z_, g_] :=
  Module[{prec}, prec = precRound[Precision[{n, m, z, g}]];
     ((-1)^(m + 1)/(g*numericJF[n, -m, g, prec]*normfactor[n, m, g^2, prec]*
        normfactor[n, -m, g^2, prec]))*radialPrefactor[m,z]*numericQS[n, m, 2, g, prec, z]] /; 
   numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity && !outsideRadius[z,g] && 
    intQ[m, n] && m >= 0
SpheroidalS2Prime[n_, m_, z_, g_] := 
  Module[{prec}, prec = precRound[Precision[{n, m, z, g}]];
     ((-1)^(m + 1)/(g*numericJF[n, -m, g, prec]*normfactor[n, m, g^2, prec]*
        normfactor[n, -m, g^2, prec]))*radialPrefactor[m,z]*numericQSPrime[n, m, 2, g, prec,
      z]] /; numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && m >= 0

(* Bug workaround *)

SpheroidalS2[n_, m_, z_, g_] := SpheroidalS2[n, -m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && m < 0

SpheroidalS2Prime[n_, m_, z_, g_] := SpheroidalS2Prime[n, -m, z, g] /;
  numQ[n, m, g, z] && Precision[{n, m, g, z}] < Infinity &&
    !outsideRadius[z,g] && intQ[m, n] && m < 0

numericJF[nu_, mu_, g_, prec_] := (numericJF[nu, mu, g, prec] =
Module[{M, prec1, nu1, mu1, g1, list, A, gamma, sum1, sum2},
    M = Max[$seriesLength, Round[Abs[g]]];
    prec1 = precRound[prec + $precIncInitSmall]; While[True, 
      {mu1, nu1, g1} = SetPrecision[{mu, nu, g}, prec1]; 
       list = angcoeff[nu, -mu, g^2, prec1, M];
       A = Table[(-1)^k, {k, -M, M}] . list;
       If[A==0, prec1 += $precInc; Continue[]];
       list = Reverse[Take[list, M + 1]]; 
       gamma = Table[(-1)^k/(Gamma[-k + nu1 + 3/2]*k!), {k, 0, M}]; 
       sum1 = list . gamma/A; If[Precision[sum1] < prec, prec1 += $precInc]; 
       If[Abs[Last[list]] > Abs[sum1]/10^prec, M += $seriesInc]; 
       If[Precision[sum1] >= prec && Abs[Last[list]] < Abs[sum1]/10^prec, 
        Break[]]]; list = angcoeff[nu, mu, g^2, prec1, M]; 
     gamma = Table[(-1)^k/(Gamma[-k - nu1 + 1/2]*k!), {k, -M, M}]; 
     sum2 = list . gamma; N[Exp[nu1*Pi*I]*2^(-2*nu1 - 1)*Gamma[nu1 - mu1 + 1]*
       g1^nu1*(sum1/sum2), prec]])
SpheroidalJoiningFactor[nu_, mu_, g_] := 
  numericJF[nu, mu, g, Max[Precision[{nu, mu, g}], 16]] /; 
   Precision[{nu, mu, g}] < Infinity && notHalfInt[nu]

SpheroidalClearCache[] := (
        clearValues[Spheroidal`Private`numericJF];
        clearValues[Spheroidal`Private`bouwkamp];
        clearValues[Spheroidal`Private`angcoeff];
        clearValues[Spheroidal`Private`angcoeff2];
        clearValues[Spheroidal`Private`normfactor]; Null)

Derivative[0, 0, 0, 1, 0][SpheroidalPS][nu_, mu_, type_:1, z_, g_] := 
  SpheroidalPSPrime[nu, mu, type, z, g]
Derivative[0, 0, 0, 1, 0][SpheroidalQS][nu_, mu_, type_:1, z_, g_] := 
  SpheroidalQSPrime[nu, mu, type, z, g]
Derivative[0, 0, 0, 1, 0][SpheroidalS1][a_, nu_, mu_, z_, g_] := 
  SpheroidalS1Prime[a, nu, mu, z, g]
Derivative[0, 0, 0, 1, 0][SpheroidalS2][a_, nu_, mu_, z_, g_] :=
  SpheroidalS2Prime[a, nu, mu, z, g]
Derivative[0, 0, 0, 1, 0][SpheroidalS3][a_, nu_, mu_, z_, g_] := 
  SpheroidalS3Prime[a, nu, mu, z, g]
Derivative[0, 0, 0, r_Integer, 0][SpheroidalPS][nu_, mu_, type_:1, z_, g_] := 
  Derivative[0, 0, 0, r - 1, 0][SpheroidalPSPrime][nu, mu, type, z, g] /; 
   r > 1
Derivative[0, 0, 0, r_Integer, 0][SpheroidalQS][nu_, mu_, type_:1, z_, g_] := 
  Derivative[0, 0, 0, r - 1, 0][SpheroidalQSPrime][nu, mu, type, z, g] /; 
   r > 1
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS1][a_, nu_, mu_, z_, g_] := 
  Derivative[0, 0, 0, r - 1, 0][SpheroidalS1Prime][a, nu, mu, z, g] /; r > 1
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS2][a_, nu_, mu_, z_, g_] :=
  Derivative[0, 0, 0, r - 1, 0][SpheroidalS2Prime][a, nu, mu, z, g] /; r > 1
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS3][a_, nu_, mu_, z_, g_] := 
  Derivative[0, 0, 0, r - 1, 0][SpheroidalS3Prime][a, nu, mu, z, g] /; r > 1
Derivative[0, 0, 0, r_Integer, 0][SpheroidalPSPrime][nu_, mu_, type_:1, z_, 
   g_] := D[(1/(1 - z^2))*(2*z*SpheroidalPSPrime[nu, mu, type, z, g] - 
     (SpheroidalEigenvalue[nu, mu, g] + g^2*(1 - z^2) - mu^2/(1 - z^2))*
      SpheroidalPS[nu, mu, type, z, g]), {z, r - 1}]
Derivative[0, 0, 0, r_Integer, 0][SpheroidalQSPrime][nu_, mu_, type_:1, z_, 
   g_] := D[(1/(1 - z^2))*(2*z*SpheroidalQSPrime[nu, mu, type, z, g] - 
     (SpheroidalEigenvalue[nu, mu, g] + g^2*(1 - z^2) - mu^2/(1 - z^2))*
      SpheroidalQS[nu, mu, type, z, g]), {z, r - 1}]
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS1Prime][a_, nu_, mu_, z_, g_] := 
  D[(1/(1 - z^2))*(2*z*SpheroidalS1Prime[a, nu, mu, z, g] - 
     (SpheroidalEigenvalue[nu, mu, g] + g^2*(1 - z^2) - mu^2/(1 - z^2))*
      SpheroidalS1[a, nu, mu, z, g]), {z, r - 1}]
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS2Prime][a_, nu_, mu_, z_, g_] := 
  D[(1/(1 - z^2))*(2*z*SpheroidalS2Prime[a, nu, mu, z, g] - 
     (SpheroidalEigenvalue[nu, mu, g] + g^2*(1 - z^2) - mu^2/(1 - z^2))*
      SpheroidalS2[a, nu, mu, z, g]), {z, r - 1}]
Derivative[0, 0, 0, r_Integer, 0][SpheroidalS3Prime][a_, nu_, mu_, z_, g_] := 
  D[(1/(1 - z^2))*(2*z*SpheroidalS3Prime[a, nu, mu, z, g] - 
     (SpheroidalEigenvalue[nu, mu, g] + g^2*(1 - z^2) - mu^2/(1 - z^2))*
      SpheroidalS3[a, nu, mu, z, g]), {z, r - 1}]

psA[nu_, mu_, k_] := recA[nu, mu, k/2, 1]
psB[nu_, mu_, k_] := -2*(((nu + k)*(nu + k + 1) + mu^2 - 1)/
    ((2*nu + 2*k - 1)*(2*nu + 2*k + 3)))
psC[nu_, mu_, k_] := recC[nu, mu, k/2, 1]
psa[nu_, mu_, j_Integer, 0] = 0; 
psa[nu_, mu_, 0, 0] = 1; 
psa[nu_, mu_, j_, k_] := 0 /; j < Abs[k]/2
pslam[nu_, mu_, (j_)?Negative] := 0
pslam[nu_, mu_, 0] = nu*(nu + 1); 
pslam[nu_, mu_, 1] = -2*((nu*(nu + 1) + mu^2 - 1)/((2*nu - 1)*(2*nu + 3))); 
pslam[nu_, mu_, j_] := psa[nu, mu, j - 1, -2]*psA[nu, mu, 0] + psa[nu, mu, j - 1, 2]*psC[nu, mu, 0]
psa[nu_, mu_, j_, k_] := psa[nu,mu,j,k] = (1/(k*(2*nu + k + 1)))*
    (Sum[pslam[nu, mu, i]*psa[nu, mu, j - i, k], {i, 1, j - Abs[k]/2}] - 
     psa[nu, mu, j - 1, k - 2]*psA[nu, mu, k] - psa[nu, mu, j - 1, k]*
      psB[nu, mu, k] - psa[nu, mu, j - 1, k + 2]*psC[nu, mu, k])

pspoch[nu_, mu_, k_] := Module[{p1 = Pochhammer[nu + mu + 1, k], p2 = Pochhammer[nu - mu + 1, k]},
                        If[p1 == 0 || p2 == 0 ||
                           p1 === ComplexInfinity || p2 === ComplexInfinity, 0, (2*nu + 1)*(p1/((2*nu + 2*k + 1)*p2))]]
psa0[nu_, mu_, t_, g_] := 
  (Sum[Sum[psa[nu, mu, j, 2*k]*g^(2*j), {j, Abs[k], t}]^2*
      pspoch[nu, mu, 2*k], {k, -t, t}] + O[g]^(2*t + 2))^(-2^(-1))

psC2[nu_, mu_] := Module[{del = Mod[mu + nu, 2]}, (-1)^del/((2*mu - 2*del - 1)*(2*mu - 2*del + 1))]
psa2[nu_, mu_, j_, k_] := 0 /; j < Abs[k]/2
psa2[nu_, mu_, j_, k_] := (1/(k*(2*nu + k + 1)))*(Sum[pslam[nu, mu, i]*psa2[nu, mu, j - i, k], 
       {i, 1, j - Abs[k]/2}] - psa2[nu, mu, j - 1, k - 2]*psA[nu, mu, k] - 
      psa2[nu, mu, j - 1, k]*psB[nu, mu, k] - psa2[nu, mu, j - 1, k + 2]*
       psC[nu, mu, k]) /; k < -mu - nu - 2 + Mod[mu + nu, 2]
psa2[nu_, mu_, j_, k_] := (1/(k*(2*nu + k + 1)))*(Sum[pslam[nu, mu, i]*psa2[nu, mu, j - i, k], 
       {i, 1, j - Abs[k]/2}] - psa2[nu, mu, j - 1, k - 2]*psA[nu, mu, k] - 
      psa2[nu, mu, j - 1, k]*psB[nu, mu, k] - psa[nu, mu, j - 1, k + 2]*
       psC2[nu, mu]) /; k == -mu - nu - 2 + Mod[mu + nu, 2]

NaturalQ = Function[ord, IntegerQ[ord] && ord > 0]

Derivative[0, 0, ord_?NaturalQ][SpheroidalEigenvalue][nu_, mu_, 0] := If[OddQ[ord], 0, ord! pslam[nu, mu, ord/2]]

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalPS][nu_, mu_, type_:1, z_, 0] :=  If[OddQ[ord], 0,
   FullSimplify[Module[{g, t=ord/2}, D[psa0[nu, mu, t, g]*Sum[Sum[psa[nu, mu, j, 2*k]*LegendreP[nu + 2*k, mu, type, z]*(-1)^k,
               {k, -j, j}]*g^(2*j), {j, 0, t}], {g,ord}]/.{g->0}]]] /; typederiv == 0

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalQS][n_, m_, type_:1, z_, 0] :=  If[OddQ[ord], 0,
   FullSimplify[Module[{g, t=ord/2}, D[psa0[n, m, t, g]*Sum[Sum[If[2*k > -Abs[m] - n - 1,
          psa[n, m, j, 2*k]*LegendreQ[n + 2*k, m, type+1, z]*(-1)^k,
          psa2[n, m, j, 2*k]*LegendreP[-n - 2*k - 1, m, type+1, z]*(-1)^k],
          {k, -j, j}]*g^(2*j), {j, 0, t}], {g,ord}]/.{g->0}]]] /; intQ[m, n] && typederiv == 0

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalQS][nu_, mu_, type_:1, z_, 0] :=  If[OddQ[ord], 0,
   FullSimplify[Module[{g, t=ord/2}, D[psa0[nu, mu, t, g]*Sum[Sum[psa[nu, mu, j, 2*k]*LegendreQ[nu + 2*k, mu, type+1, z]*(-1)^k,
          {k, -j, j}]*g^(2*j), {j, 0, t}], {g,ord}]/.{g->0}]]] /; !intQ[mu + nu]

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalQS][nu_, mu_, type_:1, z_, 0] :=  If[OddQ[ord], 0,
   FullSimplify[Module[{g, t=ord/2, expterm = If[type == 1, Cos[mu*Pi], Exp[I*mu*Pi]]},
         D[psa0[nu, mu, t, g]*Sum[Sum[If[2*k > -mu - nu - 1, 
          psa[nu, mu, j, 2*k]*LegendreQ[nu + 2*k, mu, type+1, z]*(-1)^k,
          psa2[nu, mu, j, 2*k]*(expterm*Cos[Pi*nu]*LegendreP[-nu - 2*k - 1, 
              mu, type+1, z] + (Sin[Pi*(nu - mu)]/Pi)*LegendreQ[
              -nu - 2*k - 1, mu, type+1, z])], {k, -j, j}]*g^(2*j), {j, 0, t}], {g,ord}]/.{g->0}]]] /; intQ[mu + nu] && !intQ[mu, nu] && typederiv == 0

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalPSPrime][nu_, mu_, type_:1, z_, 0] := Module[{zz}, FullSimplify[D[Derivative[0, 0, 0, ord][SpheroidalPS][nu, mu, type, zz, 0], zz]/.{zz->z}]] /; typederiv == 0
Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalQSPrime][nu_, mu_, type_:1, z_, 0] := Module[{zz}, FullSimplify[D[Derivative[0, 0, 0, ord][SpheroidalQS][nu, mu, type, zz, 0], zz]/.{zz->z}]] /; typederiv == 0

SpheroidalS1[nu_, mu_, type_:1, z_, g_] := KroneckerDelta[nu] /; g==0
SpheroidalS1Prime[nu_, mu_, type_:1, z_, g_] := 0 /; g==0
                                 
TopS1[nu_, mu_,z_,ord_?NaturalQ] := Module[{t=ord/2}, Sum[Sum[psa[nu, -Abs[mu], j, 2*k]*Series[besJ[(nu+2*k), g*z],{g,0,ord}], {k, -j, j}]
                                                              *g^(2*j), {j, 0, t}]*(1-1/z^2)^(Abs[mu]/2)]
                                 
BottomS1[nu_, mu_,z_,ord_?NaturalQ] := Module[{t=ord/2}, Sum[ Sum[psa[nu, -Abs[mu], j, 2*k]*(-1)^k,{k, -j, j}]*g^(2*j), {j, 0, t}]]                                
OrdList[ord_?NaturalQ] := Array[1+# &, ord]

U1[nu_, mu_,z_,ord_?NaturalQ] := If[EvenQ[nu],CoefficientList[TopS1[nu, mu,z,ord],g],CoefficientList[TopS1[nu, mu,z,ord],g][[OrdList[ord]]]]   
V[nu_, mu_,z_,ord_?NaturalQ] := Reverse[CoefficientList[Series[BottomS1[nu, mu,z,ord]^-1, {g, 0, ord}],g]]
U[nu_, mu_,z_,ord_?NaturalQ] := If[Length[U1[nu, mu,z,ord]]==Length[V[nu, mu,z,ord]],U1[nu, mu,z,ord],Array[0*# &, Length[V[nu, mu,z,ord]]]]                                    
                                     
Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalS1][nu_, mu_, type_:1, zz_, 0] := If[OddQ[ord+nu], 0,
   Simplify[U[nu,mu,z,ord].V[nu, mu,z,ord]*(ord!)]/.{z->zz}]  /; typederiv == 0 && intQ[nu, mu]

Derivative[0, 0, typederiv_:0, 0, ord_?NaturalQ][SpheroidalS1Prime][nu_, mu_, type_:1, z_, 0] := Module[{zz}, FullSimplify[D[Derivative[0, 0, 0, ord][SpheroidalS1][nu, mu, type, zz, 0], zz]/.{zz->z}]]  /; typederiv == 0

SpheroidalS2 /: Series[SpheroidalS2[nu_, m_, type_:1, z_, g_], {g_, 0, o_?IntegerQ}] :=
    Simplify[Module[{zz, t, mu = Abs[m], ord = o+nu+If[OddQ[o+nu], 1, 2]}, t=ord/2;
                      (2^(2*nu+1)*Exp[-(mu+nu+1)*I*Pi]*g^(-nu-1)/Gamma[nu+mu+1]*
                      Series[SpheroidalQS[nu, mu, type+1, zz, g], {g, 0, ord}]*
                      (Sum[Sum[psa[nu, -mu, j, 2*k]*(-1)^k/(Gamma[-k-nu+1/2]*k!), {k, 0, j}]*g^(2*j), {j, 0, t}]/
                       Sum[Sum[psa[nu, -mu, j, 2*k]*(-1)^k, {k, -j, j}]*g^(2*j), {j, 0, t}])/
                      (psa0[nu, mu, t, g]*Sum[Sum[psa[nu, mu, j, -2*k]*(-1)^k/(Gamma[-k+nu+3/2]*k!), {k,  0, j}]*g^(2*j), {j, 0, t}])
                      )/.{zz->z}]+ O[g]^(o+1)] /; intQ[nu, m]

SpheroidalS2 /: Series[SpheroidalS2[nu_, m_, type_:1, z_, g_], {gg_, 0, o_?IntegerQ}] :=
    (Normal[Series[SpheroidalS2[nu,m,type,z,ggg],{ggg,0,o}]]/.{ggg->g}) + O[gg]^(o+1) /; !FreeQ[g,gg]

SpheroidalS3 /: Series[SpheroidalS3[nu_, mu_, type_:1, z_, g_], {gg_, 0, o_?IntegerQ}] := 
    Simplify[Series[SpheroidalS1[nu, mu, type, z, g], {gg, 0, o}]
                 + I Series[SpheroidalS2[nu, mu, type, z, g], {gg, 0, o}]]

SpheroidalS2Prime /: Series[SpheroidalS2Prime[nu_, mu_, type_:1, z_, g_], {gg_, 0, ord_}] :=
    Simplify[Module[{zz}, D[Series[SpheroidalS2[nu, mu, type, zz, g], {gg, 0, ord}], zz] /. zz -> z]]

SpheroidalS3Prime /: Series[SpheroidalS3Prime[nu_, mu_, type_:1, z_, g_], {gg_, 0, ord_}] :=
    Simplify[Module[{zz}, D[Series[SpheroidalS3[nu, mu, type, zz, g], {gg, 0, ord}], zz] /. zz -> z]]

aspA[l_, m_, k_] = -4^(-1); 
aspB[l_, m_, k_] = m; 
aspC[l_, m_, k_] = -(2*l + 2*k + 1); 
aspD[l_, m_, k_] = (2*(l + k)^2 + 2*(l + k) - 4*m^2 + 3)/4; 
aspE[l_, m_, k_] = (-m)*(l + k + 1)*(l + k + 2); 
aspF[l_, m_, k_] = (-(l + k + 1))*(l + k + 2)*(l + k + 3)*((l + k + 4)/4); 
aspb[l_, m_, 0, 0] = 1; 
aspb[l_, m_, j_Integer, 0] = 0; 
aspb[l_, m_, j_, k_] := 0 /; Abs[k] > 4*j
aspb[l_, m_, j_, k_] := aspb[l, m, j, k] = 
   (-(aspC[l, m, k] - aspC[l, m, 0])^(-1))*
    (aspb[l, m, j - 1, k - 4]*aspA[l, m, k] + aspb[l, m, j - 1, k - 2]*
      aspB[l, m, k] + aspb[l, m, j - 1, k]*(aspD[l, m, k] - aspD[l, m, 0]) + 
     aspb[l, m, j - 1, k + 2]*aspE[l, m, k] + aspb[l, m, j - 1, k + 4]*
      aspF[l, m, k] + Sum[aspg[l, m, j - i - 1]*aspb[l, m, i, k], 
      {i, 1, j - 2}])
aspg[l_, m_, j_?Negative] := If[j==-2,-1,If[j==-1,-aspC[l, m, 0]],0]
aspg[l_, m_, 0] = -aspD[l, m, 0]; 
aspg[l_, m_, j_] := aspg[l, m, j] = -(aspA[l, m, 0]*aspb[l, m, j, -4] + 
     aspB[l, m, 0]*aspb[l, m, j, -2] + aspE[l, m, 0]*aspb[l, m, j, 2] + 
     aspF[l, m, 0]*aspb[l, m, j, 4])
SpheroidalEigenvalue /: Series[SpheroidalEigenvalue[n_, m_, g_], 
    {g_, Infinity, ord_Integer}] := 
   Module[{series, len = Length[DownValues[aspb]]}, 
     series = If[ord < -2, SeriesData[g, Infinity, {}, 0, 2, 1] /. g -> 1/g, 
        SeriesData[g, 0, Table[aspg[n - m, m, j], {j, -2, ord + 1}], -2, 
          ord + 1, 1] /. g -> 1/g]; clearValue[aspb, len]; series] /; 
    If[numQ[m, n], intQ[m, n] && n >= m >= 0, True]
SpheroidalEigenvalue /: Series[SpheroidalEigenvalue[n_, m_, g_], 
    {g_, -Infinity, ord_Integer}] := 
   Module[{series, len = Length[DownValues[aspb]]}, 
     series = If[ord < -2, SeriesData[g, Infinity, {}, 0, 2, 1] /. g -> 1/g, 
        SeriesData[g, 0, Table[(-1)^j*aspg[n - m, m, j], {j, -2, ord + 1}], 
          -2, ord + 1, 1] /. g -> 1/g]; clearValue[aspb, len]; series] /; 
    If[numQ[m, n], intQ[m, n] && n >= m >= 0, True]

asoA[s_, m_, k_] = (s + k + 1)*(m + s + k + 1); 
asoB[s_, m_, k_] = -2*k*(m + 2*s + k + 1); 
asoC[s_, m_, k_] = (s + k)*(m + s + k); 
asos[s_, m_, 0, 0] = 1; 
asos[s_, m_, j_Integer, 0] = 0; 
asos[s_, m_, j_, k_] := 0 /; Abs[k] > j
asos[s_, m_, j_, k_] := (Sum[asoev[s, m, i]*asos[s, m, j - i - 1, k], {i, 1, j - 2}] - 
     asoA[s, m, k]*asos[s, m, j - 1, k + 1] - asoB[s, m, k]*
      asos[s, m, j - 1, k] - asoC[s, m, k]*asos[s, m, j - 1, k - 1])/(4*k)
asoev[s_, m_, (j_)?Negative] = 0; 
asoev[s_, m_, -1] = 2*(m + 2*s + 1); 
asoev[s_, m_, 0] = -2*s*(m + s + 1) - m - 1; 
asoev[s_, m_, j_] := asoA[s, m, 0]*asos[s, m, j, 1] + asoC[s, m, 0]*asos[s, m, j, -1]
SpheroidalEigenvalue /: Series[SpheroidalEigenvalue[n_, m_, g_], 
    {g_, (-I)*Infinity, ord_Integer}] := 
   Module[{s = If[OddQ[n - m], n - m - 1, n - m]/2, 
      len = Length[DownValues[asos]]}, 
     series = If[ord < -1, SeriesData[g, Infinity, {}, 0, 1, 1] /. g -> 1/g, 
        SeriesData[g, 0, Table[(-I)^j*asoev[s, m, j], {j, -1, ord + 1}], -1, 
          ord + 1, 1] /. g -> 1/g]; clearValue[asos, len]; series] /; 
    If[numQ[m, n], intQ[m, n] && n >= m >= 0, True]
SpheroidalEigenvalue /: Series[SpheroidalEigenvalue[n_, m_, g_], 
    {g_, I*Infinity, ord_Integer}] := 
   Module[{s = If[OddQ[n - m], n - m - 1, n - m]/2, 
      len = Length[DownValues[asos]]}, 
     series = If[ord < -2, SeriesData[g, Infinity, {}, 0, 1, 1] /. g -> 1/g, 
        SeriesData[g, 0, Table[I^j*asoev[s, m, j], {j, -1, ord + 1}], -1, 
          ord + 1, 1] /. g -> 1/g]; clearValue[asos, len]; series] /;
    If[numQ[m, n], intQ[m, n] && n >= m >= 0, True]

boostPrecision[x_, minprec_] := SetPrecision[x, Max[Precision[x], minprec]]
    
s3[n_, m_, z_, c_, s1_, s2_] :=
  Block[{n1, m1, z1, c1, count=0,
    $seriesLength = Max[$seriesLength, Round[Abs[c]]],
    res = s1[n, m, z, c] + I s2[n, m, z, c],
    prec = precRound[Precision[{n, m, z, c}]], prec1=prec},
While[True,
    If[Precision[res] >= prec, Break[],
    prec1 += $precInc;
    count++;
    n1 = boostPrecision[n, prec1];
    m1 = boostPrecision[m, prec1];
    z1 = boostPrecision[z, prec1];
    c1 = boostPrecision[c, prec1];
    Block[{$seriesLength = $seriesLength + count*$seriesInc},
    res = s1[n1, m1, z1, c1] + I s2[n1, m1, z1, c1]]]]; res]
    
SpheroidalS3[n_, m_, z_, c_] := s3[n, m, z, c, SpheroidalS1, SpheroidalS2] /; numQ[n,m,z,c]
SpheroidalS3Prime[n_, m_, z_, c_] := s3[n, m, z, c, SpheroidalS1Prime, SpheroidalS2Prime] /; numQ[n,m,z,c]

friendlyOn[General::spell]; 
friendlyOn[General::spell1]; 

End[]; 

SetAttributes[{SpheroidalEigenvalue, SpheroidalPS, SpheroidalPSPrime, 
    SpheroidalQS, SpheroidalQSPrime, SpheroidalRadialFactor,
    SpheroidalJoiningFactor, SpheroidalS1, SpheroidalS1Prime, SpheroidalS2,
    SpheroidalS2Prime, SpheroidalS3, SpheroidalS3Prime}, 
     {Listable, Protected, ReadProtected, NumericFunction}];

EndPackage[];
