BeginPackage["DisturbingFunction`"]

(*--------------------------Version 1.2---------------------------------*) 
 (*-------------Carl Murray,  QMW,  London,  
      October 1999---------- - ----*)

DisturbingFunction::usage "DisturbingFunction.m is a package that \
defines a lot of useful functions for use in the derivation of the planetary \
disturbing function."

\(ValidArguments::usage = \*""ValidArguments[c1,c2,pw] produces an ordered list of valid arguments of the expansion of the disturbing function up to and including degree 'pw' which start with (Subscript[c, 1])λ'+(Subscript[c, 2])λ.  The components of each argument, {(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[c, 5]),(Subscript[c, 6])}, represent the argument (Subscript[c, 1])λ'+(Subscript[c, 2])λ+(Subscript[c, 3])ϖ'+(Subscript[c, 4])ϖ+(Subscript[c, 5])Ω'+(Subscript[c, 6])Ω.  c1 and c2 are integers of the form, e.g. j and i-j, where i, the order of the resonance, is an integer.""\)

\(ValidArgumentsE::usage = \*""ValidArgumentsE[(Subscript[c, 1]),(Subscript[c, 2]),pw] produces an ordered list of valid arguments of the expansion of the indirect part (external perturber) of the disturbing function up to and including degree 'pw' which start with (Subscript[c, 1])λ'+(Subscript[c, 2])λ.  The components of each argument, {(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[c, 5]),(Subscript[c, 6])}, represent the argument (Subscript[c, 1])λ'+(Subscript[c, 2])λ+(Subscript[c, 3])ϖ'+(Subscript[c, 4])ϖ+(Subscript[c, 5])Ω'+(Subscript[c, 6])Ω.  (Subscript[c, 1]) and (Subscript[c, 2]) are integers of the form, e.g. j and i-j, where i, the order of the resonance, is an integer.""\)

\(ValidArgumentsI::usage = \*""ValidArgumentsI[(Subscript[c, 1]),(Subscript[c, 2]),pw] produces an ordered list of valid arguments of the expansion of the indirect part(internal perturber) of the disturbing function up to and including degree 'pw' which start with (Subscript[c, 1])λ'+(Subscript[c, 2])λ.  The components of each argument, {(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[c, 5]),(Subscript[c, 6])}, represent the argument (Subscript[c, 1])λ'+(Subscript[c, 2])λ+(Subscript[c, 3])ϖ'+(Subscript[c, 4])ϖ+(Subscript[c, 5])Ω'+(Subscript[c, 6])Ω.  (Subscript[c, 1]) and (Subscript[c, 2]) are integers of the form, e.g. j and i-j, where i, the order of the resonance, is an integer.""\)

\(IndirectTestE::usage = \*""IndirectTestE[{(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[c, 5]),(Subscript[c, 6])}] tests whether or not the argument (Subscript[c, 1])λ'+(Subscript[c, 2])λ+(Subscript[c, 3])ϖ'+(Subscript[c, 4])ϖ+(Subscript[c, 5])Ω'+(Subscript[c, 6])Ω is a valid argument of the expansion of the indirect part (external perturber) of the expansion of the disturbing function.  All the (Subscript[c, i])'s must be explicit integers.""\)

\(IndirectTestI::usage = \*""IndirectTestI[{(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[c, 5]),(Subscript[c, 6])}] tests whether or not the argument (Subscript[c, 1])λ'+(Subscript[c, 2])λ+(Subscript[c, 3])ϖ'+(Subscript[c, 4])ϖ+(Subscript[c, 5])Ω'+(Subscript[c, 6])Ω is a valid argument of the expansion of the indirect part (internal perturber) of the expansion of the disturbing function.  All the (Subscript[c, i])'s must be explicit integers.""\)

FlmpAllan::usage = "FlmpAllan[l,m,p,inc,pw] gives the series \
expansion of Allan's F(l,m,p) inclination function up to and including powers \
of pw in inc."

FlmpKaula::usage = "FlmpKaula[l,m,p,inc,pw] gives the series \
expansion of Kaula's F(l,m,p) inclination function up to and including powers \
of pw in inc."

Newcomb::usage = "Newcomb[a,b,c,d] returns the Newcomb operator, \
usually denoted by X(a,b,c,d) which is a coefficient in the Hansen series \
expansion in the eccentricity."

Hansen::usage = "Hansen[a,b,c,ecc,pw] gives the series for the Hansen \
coefficient X(a,b,c) eccentricity function up to powers of pw in ecc.  The \
coefficient of each term in the series is the Newcomb operator, \
X(a,b,c,d)."

LaplaceSeries::usage = "LaplaceSeries[s,j,a,n,pw] gives the nth \
derivative of the Laplace coefficient b(s,j)(a) as a series in a up to and \
including terms of order pw.  s is a half-integer (e.g. 1/2, 3/2); j, n and \
pw are integers."

LaplaceCoefficient::usage = "LaplaceCoefficient[s,j,a,n] gives the \
nth derivative of the Laplace coefficient b(s,j)(a) as a numerical value. s \
is a half-integer (e.g. 1/2, 3/2); j, n are integers."

SMultiply::usage = "SMultiply[s1,s2,max] multiplies s1 by s2 \
truncating powers of e, e', s, s' to order max."

EllisKaulaD::usage = "EllisKaulaD[argument,pmax] finds the main \
coefficient of the direct part of the disturbing function which contains the \
cosine of argument explicitly."

EllisKaulaE::usage = "EllisKaulaE[argument,pmax] finds the main \
coefficient of the indirect part of the disturbing function (for an external \
perturber) which contains the cosine of argument explicitly."

EllisKaulaI::usage = "EllisKaulaI[argument,pmax] finds the main \
coefficient of the indirect part of the disturbing function (for an internal \
perturber) which contains the cosine of argument explicitly."

DirectPart::usage = "DirectPart[argument,pmax] finds the main \
coefficient of the direct part of the disturbing function which contains the \
cosine of argument, or its negative."

IndirectPartE::usage = "IndirectPartE[argument,pmax] finds the main \
coefficient of the indirect part of the disturbing function which contains \
the cosine of argument, or its negative, for an external perturber."

IndirectPartI::usage = "IndirectPartI[argument,pmax] finds the main \
coefficient of the indirect part of the disturbing function which contains \
the cosine of argument, or its negative, for an internal perturber."

DisturbingFunctionE::usage = "DisturbingFunctionE[argument,pmax] \
finds the full coefficient of the disturbing function (direct and indirect \
parts) which contains the cosine of argument for an external \
perturber."

DisturbingFunctionI::usage = "DisturbingFunctionI[argument,pmax] \
finds the full coefficient of the disturbing function (direct and indirect \
parts) which contains the cosine of argument for an internal \
perturber."

\(BigRE::usage = \*""BigRE[(Subscript[c, 1]),(Subscript[c, 2]),pw] finds the expansion of the disturbing function for an external perturber when the basic form of the argument is (Subscript[c, 1])λ'+(Subscript[c, 2])λ+ ... .  The result will include the secular terms (i.e. those independent of lamda' and lamda), and those resonant terms associated with (Subscript[c, 1])λ'+(Subscript[c, 2])λ + ... and integer multiples of the argument up to and including terms which are of degree ≤ pw in the eccentricities and inclinations.""\)

\(NBigRE::usage = \*""NBigRE[(Subscript[c, 1]),(Subscript[c, 2]),pw,α] finds the expansion of the disturbing
	function for an external perturber when the basic form of the
	argument is (Subscript[c, 1])λ'+(Subscript[c, 2])λ+ ... and calculates the numerical
	values of the Laplace coefficients taking a/a' = α.""\)

\(BigRI::usage = \*""BigRI[(Subscript[c, 1]),(Subscript[c, 2]),pw] finds the expansion of the disturbing function for an internal perturber when the basic form of the argument is (Subscript[c, 1])λ'+(Subscript[c, 2])λ+ ... .  The result will include the secular terms (i.e. those independent of lamda' and lamda), and those resonant terms associated with (Subscript[c, 1])λ'+(Subscript[c, 2])λ+ ... and integer multiples of the argument up to and including terms which are of degree ≤ pw in the eccentricities and inclinations.""\)

\(NBigRI::usage = \*""NBigRI[(Subscript[c, 1]),(Subscript[c, 2]),pw,α] finds the expansion of the disturbing
	function for an internal perturber when the basic form of the
	argument is (Subscript[c, 1])λ'+(Subscript[c, 2])λ+ ... and calculates the numerical
	values of the Laplace coefficients taking a/a' = α.""\)

Global`b /: Global`b[sss_, jjj_, nd_] := 
        Global`b[sss, -jjj, nd] /; jjj < 0;

Begin["`Private`"]

ValidArguments[c1_, c2_, n_Integer] := 
      Module[{absarglist, arglist, degree, temparg, order = Abs[c1 + c2], 
          ac1 = c1, ac2 = c2, ii, jj, kk, ll, var}, 
        arglist = 
          Flatten[Table[{kk, ll, ii, jj}, {ii, 
                If[order == 0, 0, -n], n}, {jj, 
                If[order == 0 && ii == 0, 0, -n + Abs[ii]], 
                n - Abs[ii]}, {kk, 
                If[order == 0 && ii == 0 && jj == 0, 
                  0, -n + Abs[ii] + Abs[jj]], 
                n - Abs[ii] - Abs[jj]}, {ll, 
                If[order == 0 && ii == 0 && jj == 0 && 
                    kk == 0, 0, -n + Abs[ii] + Abs[jj] + Abs[kk]], 
                n - Abs[ii] - Abs[jj] - Abs[kk]}], 
            3]; 
arglist = 
          Select[arglist, (order + Apply[Plus, #] == 
                  0) &]; 
arglist = 
          Select[arglist, 
            EvenQ[Abs[#[[1]] + #[[2]]]] &]; 
absarglist = 
          Map[Abs, arglist]; 
temparg = 
          arglist /. {x1_, x2_, x3_, x4_} :> {x3, x4, x1, 
                x2}; 
degree = 
          Apply[Plus, absarglist, 1]; 
arglist = 
          Prepend[Transpose[arglist], temparg]; 
arglist = 
          Prepend[arglist, absarglist]; 
arglist = 
          Sort[Transpose[Prepend[arglist, degree]]]; 
var = 
          Variables[{c1, c2}]; 
If[var =!= {}, 
          If[Coefficient[c1, var[[1]], 1] == -1, 
            ac1 = -c1;  ac2 = -c2], 
          If[c1 < 0, ac1 = -c1;  ac2 = -c2]]; 
If[
          Length[arglist] =!= 0, 
          Transpose[
              Drop[Transpose[arglist], 3]] /. {x1_, x2_, x3_, 
                x4_} :> {ac1, ac2, x3, x4, x1, x2}, 
          arglist]]

IndirectTestE[{clamdap_, clamda_, clpp_Integer, clp_Integer, 
          clanp_Integer, clan_Integer}] := 
      Module[{qp = -clpp, q = clp, m1pp = clanp, pq1p = clamdap - 1, 
          pq1 = -clamda - 1, p, pp, m}, 
        pp = -(pq1p - qp)/
            2; 
p = -(pq1 - q)/
            2; 
m = m1pp - 2*pp + 1; 
If[
          p ≤ 1 && p ≥ 0 && pp ≤ 1 && 
            pp ≥ 0 && m ≤ 1 && m ≥ 0 && 
            clamdap =!= 0 && IntegerQ[p] == True && 
            IntegerQ[pp] == True, True, False]]

ValidArgumentsE[c1_, c2_, pmax_Integer] := 
      Module[{valist, valen, arglist = {}, temp, inonzero, argument1, 
          argument2, var, novar}, 
        var = Variables[{c1, c2}]; 
novar = 
          If[var == {}, True, False]; 
valist = 
          ValidArguments[c1, c2, pmax]; 
valen = 
          Length[valist]; 
Do[
          temp = {}; 
If[i == 0 || novar, 
            inonzero = False, 
            inonzero = True]; 
argument1 = 
            valist[[k]]; 
If[! novar, 
            argument1 = argument1 /. var[[1]] -> i];  
          If[inonzero, 
            argument2 = -valist[[k]]]; 
If[! novar, argument2 = 
              argument2 /. 
                var[[1]] -> -i]; 
If[
            IndirectTestE[argument1] || IndirectTestE[-argument1], 
            AppendTo[temp, argument1]]; 
If[inonzero, 
            If[IndirectTestE[argument2] || IndirectTestE[-argument2], 
              AppendTo[temp, argument2]]]; 
temp = 
            Union[temp]; 
arglist = 
            Join[arglist, temp], {k, 1, valen}, {i, 1, 
            Floor[(pmax + Abs[c1 + c2])/2] + 1}]; 
If[
          novar, Union[arglist], arglist]]

IndirectTestI[{clamdap_, clamda_, clpp_Integer, clp_Integer, 
          clanp_Integer, clan_Integer}] := 
      Module[{qp = -clpp, q = clp, m1pp = clanp, pq1p = clamdap - 1, 
          pq1 = -clamda - 1, p, pp, m}, 
        pp = -(pq1p - qp)/
            2; 
p = -(pq1 - q)/
            2; 
m = m1pp - 2*pp + 1; 
If[
          p ≤ 1 && p ≥ 0 && pp ≤ 1 && 
            pp ≥ 0 && m ≤ 1 && m ≥ 0 && 
            clamda =!= 0 && IntegerQ[p] == True && 
            IntegerQ[pp] == True, True, False]]

ValidArgumentsI[c1_, c2_, pmax_Integer] := 
      Module[{valist, valen, arglist = {}, temp, inonzero, argument1, 
          argument2, var, novar}, 
        var = Variables[{c1, c2}]; 
novar = 
          If[var == {}, True, False]; 
valist = 
          ValidArguments[c1, c2, pmax]; 
valen = 
          Length[valist]; 
Do[
          temp = {}; 
If[i == 0 || novar, 
            inonzero = False, 
            inonzero = True]; 
argument1 = 
            valist[[k]]; 
If[! novar, 
            argument1 = argument1 /. var[[1]] -> i];  
          If[inonzero, 
            argument2 = -valist[[k]]]; 
If[! novar, argument2 = 
              argument2 /. 
                var[[1]] -> -i]; 
If[
            IndirectTestI[argument1] || IndirectTestI[-argument1], 
            AppendTo[temp, argument1]]; 
If[inonzero, 
            If[IndirectTestI[argument2] || IndirectTestI[-argument2], 
              AppendTo[temp, argument2]]]; 
temp = 
            Union[temp]; 
arglist = 
            Join[arglist, temp], {k, 1, valen}, {i, 0, 
            Floor[(pmax + Abs[c1 + c2])/2] + 1}]; 
If[
          novar, Union[arglist], arglist]]

FlmpAllan[l_Integer, m_Integer, p_Integer, inc_, pw_Integer] := 
      Expand[Normal[
          Series[Expand[(I^(l - m)*(l + m)!)/(2^
                        l*p!*(l - p)!)*
                Sum[(-1)^k*Binomial[2*(l - p), k]*
                    Binomial[2*p, l - m - k]*inc^(m - l + 2*p + 2*k)*
                    Normal[Series[(1 - 
                              inc^2)^((3*l - m - 2*p - 2*k)/
                              2), {inc, 0, pw}]], {k, Max[0, l - m - 2*p], 
                    Min[l - m, 2*l - 2*p]}]], {inc, 0, pw}]]]

FlmpKaula[l_Integer, m_Integer, p_Integer, inc_, pw_Integer] := 
        Expand[I^(l - m)*FlmpAllan[l, m, p, inc, pw]];

Newcomb[a_, b_, c_, d_] := Newcomb[a, -b, d, c] /; d > c

Newcomb[a_, b_, c_, d_] := 0 /; c < 0 || d < 0

Newcomb[a_, b_, c_, 
        d_] := (2*(2*b - a)*Newcomb[a, b + 1, c - 1, 0] + (b - a)*
                Newcomb[a, b + 2, c - 2, 0])/(4*c) /; 
        c > 1 && d == 0

Newcomb[a_, b_, c_, d_] := 
      Block[{j}, (-(2*(2*b + a)*
                      Newcomb[a, b - 1, c, d - 1]) - (b + a)*
                  Newcomb[a, b - 2, c, d - 2] - (c - 5*d + 4 + 4*b + a)*
                  Newcomb[a, b, c - 1, d - 1] + 
                2*(c - d + b)*
                  Sum[(-1)^j*Binomial[3/2, j]*
                      Newcomb[a, b, c - j, d - j], {j, 2, d}])/(4*d) /; 
          c ≥ d && d ≥ 2]

Newcomb[a_, b_, c_, 
        d_] := (-(2*(2*b + a)*
                    Newcomb[a, b - 1, c, d - 1]) - (b + a)*
                Newcomb[a, b - 2, c, d - 2] - (c - 5*d + 4 + 4*b + a)*
                Newcomb[a, b, c - 1, d - 1])/(4*d) /; 
        c ≥ d && d < 2 && d ≠ 0

Newcomb[a_, b_, 0, 0] := 1

Newcomb[a_, b_, 1, 0] := b - a/2

Hansen[a_, b_, c_, e_, pw_Integer] := 
      Expand[e^Abs[c - b]*
          Sum[Newcomb[a, b, sigma + Max[0, c - b], sigma + Max[0, b - c]]*
              e^(2*sigma), {sigma, 0, 
              Quotient[pw - Abs[c - b], 2]}]]

LaplaceSeries[s_, j_Integer, alpha_, n_Integer, nord_Integer] := 
      LaplaceSeries[s, -j, alpha, n, nord] /; j < 0

LaplaceSeries[s_, j_Integer, alpha_, 0, nord_Integer] := 
        Block[{ifact}, 
          If[j == 0, ifact = 2, 
            ifact = (2/j!)*Product[(s + ii), {ii, 0, j - 1}]*
                alpha^j]; 
Expand[
            ifact*(1 + 
                  Sum[Product[(s + jj - 
                              1) (s + j + jj - 
                                1)/(jj*(j + jj)), {jj, 1, ii}]*
                      alpha^(2*ii), {ii, 
                      1, (nord - j)/2}])]];

LaplaceSeries[s_, j_Integer, alpha_, n_Integer, nord_Integer] := 
      D[LaplaceSeries[s, j, alpha, 0, nord + n], {alpha, n}] /; 
        n > 0

LaplaceCoefficient[s_, jj_Integer, a_, n_Integer] := 
      Module[{outside, j = Abs[jj], cterm, cbterm, csum, i, k}, 
        outside = 
          2*Product[s + k, {k, 0, j - 1}]/j!*
            a^(j - n); 
cterm = 
          1; 
cbterm = 
          cterm*Product[j - k, {k, 0, n - 1}]; 
For[
          i = 1;  csum = 
            cbterm, (cbterm > $MachineEpsilon) || (cbterm == 
                0), i++, 
          cterm = cterm*(s + i - 1)*(s + j + i - 1)/i/(j + i)*
              a^2; 
cbterm = 
            cterm*Product[
                j + 2*i - k, {k, 0, n - 1}]; 
csum = 
            csum + cbterm]; 
Expand[
          outside*csum]]

SMultiply[s1_, s2_, max_Integer] := 
      Module[{ts1, ts2, t}, 
        ts1 = s1 /. {Global`e -> t*Global`e, 
              Global`e' -> t*Global`e', Global`s -> t*Global`s, 
              Global`s' -> t*Global`s'}; 
ts2 = 
          s2 /. {Global`e -> t*Global`e, Global`e' -> t*Global`e', 
              Global`s -> t*Global`s, 
              Global`s' -> t*Global`s'}; 
Normal[
            Series[ts1*ts2, {t, 0, max}]] /. t -> 1]

EllisKaulaD[{clamdap_, clamda_, clpp_Integer, clp_Integer, 
            clanp_Integer, clan_Integer}, pmax_Integer] := 
        Module[{qp = -clpp, q = clp, msnpp = clanp, msnp = -clan, 
            isum = 0, lmax = pmax - Abs[clan] - Abs[clanp], coeff1, coeff2, 
            coeff3, ifact, xe, xep, ifun, efun, iefun, imin, imax, nmax, 
            mmin, p, pp, fsi, fsip, lapfun, term, i, ss, n, m, abar, jj}, 
          absqqp = Abs[q] + Abs[qp]; 
If[EvenQ[msnpp], 
            mmseven = True, 
            mmseven = False]; 
ppmp = (msnpp - msnp)/
              2; 
If[ppmp < 0, pmin = -ppmp;  
            ppmin = 0, pmin = 0;  
            ppmin = ppmp]; 
sm2nmin = 
            Max[pmin, 
              ppmin, -msnp + 2*pmin, -msnpp + 
                2*ppmin]; 
ssmin = sm2nmin;  
          imin = 0; 
imax = 
            Floor[(pmax - Abs[q] - Abs[qp])/2]; 
Do[
            ifact = ((2  i)!/i!)*(((-1)^
                        i)/(2^(2  i + 1)))*(Global`α)^
                  i; 
xe = 
              Hansen[i + kk, -clamda - q, -clamda, Global`e, 
                pmax - Abs[qp] - Abs[clan] - 
                  Abs[clanp]]; 
xep = 
              Hansen[-(i + kk + 1), clamdap - qp, clamdap, Global`e', 
                pmax - Abs[q] - Abs[clan] - 
                  Abs[clanp]]; 
efun = 
              SMultiply[xe, xep, pmax]; 
Do[
              nmax = Floor[(ss - ssmin)/2]; 
Do[
                If[(EvenQ[ss] && 
                        mmseven) || (OddQ[ss] && ! mmseven), 
                  mmin = 0, mmin = 1]; 
s2mn = 
                  ss - 2*n; 
Do[
                  pp = (msnpp - m + s2mn)/
                      2; 
p = (msnp - m + s2mn)/
                      2; 
If[
                    p ≤ s2mn && p ≥ pmin && 
                      pp ≤ s2mn && pp ≥ ppmin, 
                    If[m == 0, kappa = 1, 
                      kappa = 2]; 
fsi = 
                      FlmpKaula[s2mn, m, p, Global`s, 
                        pmax - absqqp - 
                          Abs[clanp]]; 
fsip = 
                      FlmpKaula[s2mn, m, pp, Global`s', 
                        pmax - absqqp - Abs[clan]]; 
ifun = 
                      SMultiply[fsi, fsip, 
                        pmax - absqqp]; 
iefun = 
                      SMultiply[efun, ifun, 
                        pmax]; 
coeff1 = ((2*ss - 4*n + 
                                1)*(ss - n)!)/(2^(2*
                                  n)*n!*(2*ss - 2*n + 
                                  1)!); 
coeff2 = 
                      kappa*(s2mn - m)!/(s2mn + 
                                m)!; 
Do[
                      jj = clamda + i - 2*abar - 2*n - 2*p + 
                          q; 
vjj = 
                        Variables[jj]; 
If[vjj =!= {}, 
                        If[Coefficient[jj, vjj[[1]], 1] == -1, 
                          jj = -jj], 
                        jj = Abs[
                            jj]]; 
coeff3 = ((-1)^
                                ss*2^(2*ss))/((i - ss - 
                                    abar)!*abar!); 
lapfun = 
                        Sum[((-1)^l/(l!))*
                            Binomial[l, k]*(-1)^
                              k*(Global`α)^l*
                            Global`b[i + 1/2, jj, 
                              l]*(iefun /. kk -> k), {l, 0, 
                            lmax}, {k, 0, l}]; 
term = 
                        ifact*coeff1*coeff2*coeff3*
                          lapfun; 
isum = 
                        isum + term, {abar, 0, i - ss}]], {m, mmin, 
                    ss - 2  n, 2}], {n, 0, nmax}], {ss, ssmin, i}], {i, imin, 
              imax}]; 
Expand[isum]];

EllisKaulaE[{clamdap_, clamda_, clpp_Integer, clp_Integer, 
          clanp_Integer, clan_Integer}, pmax_Integer] := 
      Module[{qp = -clpp, q = clp, m1pp = clanp, m1p = -clan, 
          pq1p = clamdap - 1, pq1 = -clamda - 1, coeff, p, pp, m, kappa, 
          fsi, fsip, xe, xep, efun, ifun, iefun}, 
        pp = -(pq1p - qp)/
            2; 
p = -(pq1 - q)/
            2; 
m = m1pp - 2*pp + 1; 
If[
          p ≤ 1 && p ≥ 0 && pp ≤ 1 && 
            pp ≥ 0 && m ≤ 1 && m ≥ 0 && 
            IntegerQ[p] == True && IntegerQ[pp] == True, 
          absqqp = Abs[q] + Abs[qp]; 
If[m == 0, 
            kappa = 1, kappa = 2]; 
coeff = 
            kappa*(1 - m)!/(1 + m)!; 
fsi = 
            FlmpKaula[1, m, p, Global`s, 
              pmax - absqqp - Abs[clanp]]; 
fsip = 
            FlmpKaula[1, m, pp, Global`s', 
              pmax - absqqp - Abs[clan]]; 
ifun = 
            SMultiply[fsi, fsip, pmax - absqqp]; 
xe = 
            Hansen[1, -clamda - q, -clamda, Global`e, 
              pmax - Abs[qp] - Abs[clan] - 
                Abs[clanp]]; 
xep = 
            Hansen[-2, clamdap - qp, clamdap, Global`e', 
              pmax - Abs[q] - Abs[clan] - 
                Abs[clanp]]; 
efun = 
            SMultiply[xe, xep, pmax]; 
iefun = 
            SMultiply[efun, ifun, pmax]; 
Expand[
            coeff*iefun], 0]]

EllisKaulaI[{clamdap_, clamda_, clpp_Integer, clp_Integer, 
          clanp_Integer, clan_Integer}, pmax_Integer] := 
      Module[{qp = -clpp, q = clp, m1pp = clanp, m1p = -clan, 
          pq1p = clamdap - 1, pq1 = -clamda - 1, coeff, p, pp, m, kappa, 
          fsi, fsip, xe, xep, efun, ifun, iefun}, 
        pp = -(pq1p - qp)/
            2; 
p = -(pq1 - q)/
            2; 
m = m1pp - 2*pp + 1; 
If[
          p ≤ 1 && p ≥ 0 && pp ≤ 1 && 
            pp ≥ 0 && m ≤ 1 && m ≥ 0 && 
            IntegerQ[p] == True && IntegerQ[pp] == True, 
          absqqp = 
            Abs[q] + Abs[qp]; 
  
            If[m == 0, kappa = 1, 
            kappa = 2]; 
coeff = 
            kappa*(1 - m)!/(1 + m)!; 
fsi = 
            FlmpKaula[1, m, p, Global`s, 
              pmax - absqqp - Abs[clanp]]; 
fsip = 
            FlmpKaula[1, m, pp, Global`s', 
              pmax - absqqp - Abs[clan]]; 
ifun = 
            SMultiply[fsi, fsip, pmax - absqqp]; 
xe = 
            Hansen[-2, -clamda - q, -clamda, Global`e, 
              pmax - Abs[qp] - Abs[clan] - 
                Abs[clanp]]; 
xep = 
            Hansen[1, clamdap - qp, clamdap, Global`e', 
              pmax - Abs[q] - Abs[clan] - 
                Abs[clanp]]; 
efun = 
            SMultiply[xe, xep, pmax]; 
iefun = 
            SMultiply[efun, ifun, pmax]; 
Expand[
            coeff*iefun], 0]]

DirectPart[argument_List, pmax_Integer] := 
        If[Drop[argument, 2] =!= {0, 0, 0, 0}, 
          EllisKaulaD[argument, pmax] + EllisKaulaD[-argument, pmax], 
          If[Variables[argument] =!= {}, 
            Expand[(EllisKaulaD[argument, pmax] + 
                    EllisKaulaD[-argument, pmax])/2], 
            EllisKaulaD[argument, pmax]]];

IndirectPartE[argument_List, 
          pmax_Integer] := -EllisKaulaE[argument, pmax] - 
          EllisKaulaE[-argument, pmax];

IndirectPartI[argument_List, 
          pmax_Integer] := -EllisKaulaI[argument, pmax] - 
          EllisKaulaI[-argument, pmax];

DisturbingFunctionE[argument_List, 
          pmax_Integer] := (Global`𝒢*
              Global`m'/Global`a')*(DirectPart[argument, pmax]*
                Cos[argument[[1]]*Global`λ' + 
                    argument[[2]]*Global`λ + 
                    argument[[3]]*Global`ϖ' + 
                    argument[[4]]*Global`ϖ + 
                    argument[[5]]*Global`Ω' + 
                    argument[[6]]*Global`Ω] + 
              Global`α*IndirectPartE[argument, pmax]*
                Cos[argument[[1]]*Global`λ' + 
                    argument[[2]]*Global`λ + 
                    argument[[3]]*Global`ϖ' + 
                    argument[[4]]*Global`ϖ + 
                    argument[[5]]*Global`Ω' + 
                    argument[[6]]*
                      Global`Ω]);

DisturbingFunctionI[argument_List, 
          pmax_Integer] := (Global`𝒢*
              Global`m/Global`a')*(DirectPart[argument, pmax]*
                Cos[argument[[1]]*Global`λ' + 
                    argument[[2]]*Global`λ + 
                    argument[[3]]*Global`ϖ' + 
                    argument[[4]]*Global`ϖ + 
                    argument[[5]]*Global`Ω' + 
                    argument[[6]]*Global`Ω] + (1/
                    Global`α^2)*IndirectPartI[argument, pmax]*
                Cos[argument[[1]]*Global`λ' + 
                    argument[[2]]*Global`λ + 
                    argument[[3]]*Global`ϖ' + 
                    argument[[4]]*Global`ϖ + 
                    argument[[5]]*Global`Ω' + 
                    argument[[6]]*
                      Global`Ω]);

BigRE[c1_Integer, c2_Integer, pw_Integer] := 
        Module[{args, resargs, secargs, order = Abs[c1 + c2]}, 
          For[resargs = ValidArguments[c1, c2, pw];  i = 1, 
            i*order ≤ pw, i++, 
            resargs = 
              Union[resargs, 
                ValidArguments[i*c1, i*c2, 
                  pw]]]; 
secargs = 
            ValidArguments[0, 0, pw]; 
args = 
            Union[secargs, resargs]; 
Expand[
            Plus @@ Map[DisturbingFunctionE[#, pw] &, args]]];

NBigRE[c1_Integer, c2_Integer, pw_Integer, a_] := 
        BigRE[c1, c2, pw] /. {Global`α -> a, 
            Global`b[s_, j_, n_] -> 
              LaplaceCoefficient[s, j, a, n]};

BigRI[c1_Integer, c2_Integer, pw_Integer] := 
        Module[{args, resargs, secargs, order = Abs[c1 + c2]}, 
          For[resargs = ValidArguments[c1, c2, pw];  i = 1, 
            i*order ≤ pw, i++, 
            resargs = 
              Union[resargs, 
                ValidArguments[i*c1, i*c2, 
                  pw]]]; 
secargs = 
            ValidArguments[0, 0, pw]; 
args = 
            Union[secargs, resargs]; 
Expand[
            Plus @@ Map[DisturbingFunctionI[#, pw] &, args]]];

NBigRI[c1_Integer, c2_Integer, pw_Integer, a_] := 
        BigRI[c1, c2, pw] /. {Global`α -> a, 
            Global`b[s_, j_, n_] -> 
              LaplaceCoefficient[s, j, a, n]};

End[]

Protect[ValidArguments, IndirectTestE, IndirectTestI, ValidArgumentsE, 
      ValidArgumentsI, FlmpKaula, FlmpAllan, Newcomb, Hansen, LaplaceSeries, 
      DisturbingFunctionE, DisturbingFunctionI, SMultiply, DirectPart, 
      LaplaceCoefficient, EllisKaulaD, IndirectPartE, EllisKaulaE, 
      IndirectPartI, EllisKaulaI, BigRE, BigRI, NBigRE, NBigRI]

EndPackage[]