(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)

BeginPackage["DisturbingFunction`"]

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

(*
  Version history:
      1.0 April 1993
    1.1 March 1994 Improved version of ValidArguments
    1.2 October 1999 Improved version to use Greek letters
  *)

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 \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda].  The components of each argument, {\!\(c\_1\),\!\(c\_2\),\!\(c\_3\),\!\(c\_4\),\!\(c\_5\),\!\(c\_6\)}, represent the argument \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+\!\(c\_3\)\[CurlyPi]'+\!\(c\_4\)\[CurlyPi]+\!\(c\_5\)\[CapitalOmega]'+\!\(c\_6\)\[CapitalOmega].  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[\!\(c\_1\),\!\(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 \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda].  The components of each argument, {\!\(c\_1\),\!\(c\_2\),\!\(c\_3\),\!\(c\_4\),\!\(c\_5\),\!\(c\_6\)}, represent the argument \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+\!\(c\_3\)\[CurlyPi]'+\!\(c\_4\)\[CurlyPi]+\!\(c\_5\)\[CapitalOmega]'+\!\(c\_6\)\[CapitalOmega].  \!\(c\_1\) and \!\(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[\!\(c\_1\),\!\(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 \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda].  The components of each argument, {\!\(c\_1\),\!\(c\_2\),\!\(c\_3\),\!\(c\_4\),\!\(c\_5\),\!\(c\_6\)}, represent the argument \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+\!\(c\_3\)\[CurlyPi]'+\!\(c\_4\)\[CurlyPi]+\!\(c\_5\)\[CapitalOmega]'+\!\(c\_6\)\[CapitalOmega].  \!\(c\_1\) and \!\(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[{\!\(c\_1\),\!\(c\_2\),\!\(c\_3\),\!\(c\_4\),\!\(c\_5\),\!\(c\_6\)}] tests whether or not the argument \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+\!\(c\_3\)\[CurlyPi]'+\!\(c\_4\)\[CurlyPi]+\!\(c\_5\)\[CapitalOmega]'+\!\(c\_6\)\[CapitalOmega] is a valid argument of the expansion of the indirect part (external perturber) of the expansion of the disturbing function.  All the \!\(c\_i\)'s must be explicit integers.\>\""\)

\!\(IndirectTestI::usage = \*"\"\<IndirectTestI[{\!\(c\_1\),\!\(c\_2\),\!\(c\_3\),\!\(c\_4\),\!\(c\_5\),\!\(c\_6\)}] tests whether or not the argument \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+\!\(c\_3\)\[CurlyPi]'+\!\(c\_4\)\[CurlyPi]+\!\(c\_5\)\[CapitalOmega]'+\!\(c\_6\)\[CapitalOmega] is a valid argument of the expansion of the indirect part (internal perturber) of the expansion of the disturbing function.  All the \!\(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[\!\(c\_1\),\!\(c\_2\),pw] finds the expansion of the disturbing function for an external perturber when the basic form of the argument is \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+ ... .  The result will include the secular terms (i.e. those independent of lamda' and lamda), and those resonant terms associated with \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda] + ... and integer multiples of the argument up to and including terms which are of degree \[LessEqual] pw in the eccentricities and inclinations.\>\""\)

\!\(NBigRE::usage = \*"\"\<NBigRE[\!\(c\_1\),\!\(c\_2\),pw,\[Alpha]] finds the expansion of the disturbing\n\tfunction for an external perturber when the basic form of the\n\targument is \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+ ... and calculates the numerical\n\tvalues of the Laplace coefficients taking a/a' = \[Alpha].\>\""\)

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

\!\(NBigRI::usage = \*"\"\<NBigRI[\!\(c\_1\),\!\(c\_2\),pw,\[Alpha]] finds the expansion of the disturbing\n\tfunction for an internal perturber when the basic form of the\n\targument is \!\(c\_1\)\[Lambda]'+\!\(c\_2\)\[Lambda]+ ... and calculates the numerical\n\tvalues of the Laplace coefficients taking a/a' = \[Alpha].\>\""\)

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\[Equal]0,0,-n],n},{jj,
            If[order\[Equal]0&&ii\[Equal]0,0,-n+Abs[ii]],n-Abs[ii]},{kk,
            If[order\[Equal]0&&ii\[Equal]0&&jj\[Equal]0,0,-n+Abs[ii]+Abs[jj]],
            n-Abs[ii]-Abs[jj]},{ll,
            If[order\[Equal]0&&ii\[Equal]0&&jj\[Equal]0&&kk\[Equal]0,
              0,-n+Abs[ii]+Abs[jj]+Abs[kk]],n-Abs[ii]-Abs[jj]-Abs[kk]}],3];
    arglist=Select[arglist,(order+Apply[Plus,#]\[Equal]0)&];
    arglist=Select[arglist,EvenQ[Abs[#[[1]]+#[[2]]]]&];
    absarglist=Map[Abs,arglist];
    temparg=arglist/.{x1_,x2_,x3_,x4_}\[RuleDelayed]{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]\[Equal]-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_}\[RuleDelayed]{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\[LessEqual]1&&p\[GreaterEqual]0&&pp\[LessEqual]1&&pp\[GreaterEqual]0&&
        m\[LessEqual]1&&m\[GreaterEqual]0&&clamdap=!=0&&
        IntegerQ[p]\[Equal]True&&IntegerQ[pp]\[Equal]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\[Equal]{},True,False];
    valist=ValidArguments[c1,c2,pmax];
    valen=Length[valist];
    Do[temp={};
      If[i\[Equal]0||novar,inonzero=False,inonzero=True];
      argument1=valist[[k]];
      If[!novar,argument1=argument1/.var[[1]]\[Rule]i]; 
      If[inonzero,argument2=-valist[[k]]];
      If[!novar,argument2=argument2/.var[[1]]\[Rule]-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\[LessEqual]1&&p\[GreaterEqual]0&&pp\[LessEqual]1&&pp\[GreaterEqual]0&&
        m\[LessEqual]1&&m\[GreaterEqual]0&&clamda=!=0&&
        IntegerQ[p]\[Equal]True&&IntegerQ[pp]\[Equal]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\[Equal]{},True,False];
    valist=ValidArguments[c1,c2,pmax];
    valen=Length[valist];
    Do[temp={};
      If[i\[Equal]0||novar,inonzero=False,inonzero=True];
      argument1=valist[[k]];
      If[!novar,argument1=argument1/.var[[1]]\[Rule]i]; 
      If[inonzero,argument2=-valist[[k]]];
      If[!novar,argument2=argument2/.var[[1]]\[Rule]-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\[Equal]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\[GreaterEqual]d&&d\[GreaterEqual]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\[GreaterEqual]d&&d<2&&d\[NotEqual]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\[Equal]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\[Equal]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\[Rule]t*Global`e,Global`e'\[Rule]t*Global`e',
          Global`s\[Rule]t*Global`s,Global`s'\[Rule]t*Global`s'};
    ts2=s2/.{Global`e\[Rule]t*Global`e,Global`e'\[Rule]t*Global`e',
          Global`s\[Rule]t*Global`s,Global`s'\[Rule]t*Global`s'};
    Normal[Series[ts1*ts2,{t,0,max}]]/.t\[Rule]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=((2i)!/i!)*(((-1)^i)/(2^(2i+1)))*(Global`\[Alpha])^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\[LessEqual]s2mn&&p\[GreaterEqual]pmin&&pp\[LessEqual]s2mn&&
                  pp\[GreaterEqual]ppmin,If[m\[Equal]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]\[Equal]-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`\[Alpha])^l*
                        Global`b[i+1/2,jj,l]*(iefun/.kk\[Rule]k),{l,0,
                        lmax},{k,0,l}];
                  term=ifact*coeff1*coeff2*coeff3*lapfun;
                  isum=isum+term,{abar,0,i-ss}]],{m,mmin,ss-2n,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\[LessEqual]1&&p\[GreaterEqual]0&&pp\[LessEqual]1&&pp\[GreaterEqual]0&&
        m\[LessEqual]1&&m\[GreaterEqual]0&&IntegerQ[p]\[Equal]True&&
        IntegerQ[pp]\[Equal]True,absqqp=Abs[q]+Abs[qp];
      If[m\[Equal]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\[LessEqual]1&&p\[GreaterEqual]0&&pp\[LessEqual]1&&pp\[GreaterEqual]0&&
        m\[LessEqual]1&&m\[GreaterEqual]0&&IntegerQ[p]\[Equal]True&&
        IntegerQ[pp]\[Equal]True,absqqp=Abs[q]+Abs[qp];
      (*m=m1pp-2*pp+1;*)If[m\[Equal]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`\[ScriptCapitalG]*
          Global`m'/Global`a')*(DirectPart[argument,pmax]*
            Cos[argument[[1]]*Global`\[Lambda]'+
                argument[[2]]*Global`\[Lambda]+
                argument[[3]]*Global`\[CurlyPi]'+
                argument[[4]]*Global`\[CurlyPi]+
                argument[[5]]*Global`\[CapitalOmega]'+
                argument[[6]]*Global`\[CapitalOmega]]+
          Global`\[Alpha]*IndirectPartE[argument,pmax]*
            Cos[argument[[1]]*Global`\[Lambda]'+
                argument[[2]]*Global`\[Lambda]+
                argument[[3]]*Global`\[CurlyPi]'+
                argument[[4]]*Global`\[CurlyPi]+
                argument[[5]]*Global`\[CapitalOmega]'+
                argument[[6]]*Global`\[CapitalOmega]]);

DisturbingFunctionI[argument_List,
      pmax_Integer]:=(Global`\[ScriptCapitalG]*
          Global`m/Global`a')*(DirectPart[argument,pmax]*
            Cos[argument[[1]]*Global`\[Lambda]'+
                argument[[2]]*Global`\[Lambda]+
                argument[[3]]*Global`\[CurlyPi]'+
                argument[[4]]*Global`\[CurlyPi]+
                argument[[5]]*Global`\[CapitalOmega]'+
                argument[[6]]*Global`\[CapitalOmega]]+(1/Global`\[Alpha]^2)*
            IndirectPartI[argument,pmax]*
            Cos[argument[[1]]*Global`\[Lambda]'+
                argument[[2]]*Global`\[Lambda]+
                argument[[3]]*Global`\[CurlyPi]'+
                argument[[4]]*Global`\[CurlyPi]+
                argument[[5]]*Global`\[CapitalOmega]'+
                argument[[6]]*Global`\[CapitalOmega]]);

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\[LessEqual]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`\[Alpha]\[Rule]a,
        Global`b[s_,j_,n_]\[Rule]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\[LessEqual]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`\[Alpha]\[Rule]a,
        Global`b[s_,j_,n_]\[Rule]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[]