Chapter 6, Section 6.5

We are going to derive an expansion of the planetary disturbing function complete to second order in the the eccentricities and inclinations.

We start by defining a function, Truncate, that will truncate series in the eccentricities and inclinations to a specified order.

Truncate[expr_, n_Integer]:=
Module[{},Normal[Series[
(TrigReduce[expr] /. {e->e*p,e'->e'*p,s->s*p,s'->s'*p}),{p,0,n}]] /. p->1]

Next we use Eq.(2.84) to derive a series for Sin[f] where f is the true anomaly.

sf=
Collect[Series[TrigReduce[
2*Series[Sqrt[1-e^2],{e,0,2}]*Sum[(1/i)*
D[Normal[Series[BesselJ[i,i*e],{e,0,3}]]*
Sin[i*M],e],{i,3}]],{e,0,2}],e]

We can remove some of the brackets by collecting the terms in e.  This gives Eq.(6.73).

Collect[Expand[%],e]

Now we do a similar series for Cos[f] using Eq.(2.85).

cf=
Collect[Series[TrigReduce[
-e+(2*(1-e^2)/e)*Sum[
Normal[Series[BesselJ[i,i*e],{e,0,3}]]*
Cos[i*M],{i,3}]],{e,0,2}],e]

We can remove some of the brackets by collecting the terms in e.  This gives Eq.(6.74).

Collect[Expand[%],e]

Create similar series for Cos[f'] and Sin[f'] by substitution.

cfp=cf/.{M->M', e->e'};

sfp=sf/.{M->M', e->e'};

Now we can derive the series for Cos[ω+f], Cos[ω'+f'], Sin[ω+f] and Sin[ω'+f'].  See Eqs.(6.75) and (6.76).

cfom=Collect[TrigReduce[cf*Cos[ω]-sf*Sin[ω]],e]

cfomp=cfp*Cos[ω']-sfp*Sin[ω'];

sfom=Collect[TrigReduce[sf*Cos[ω]+cf*Sin[ω]],e]

Collect[Expand[%],e]

sfomp=sfp*Cos[ω']+cfp*Sin[ω'];

Define Cos[I] and Sin[I] in terms of s = Sin[I/2].

ci = 1 - 2  s^2;
si = 2  s;

Now we can calculate series for x/r (see Eq.(6.79)) ...

xor = 
Collect[
        Truncate[Cos[Ω]*cfom - Sin[Ω]*sfom*ci, 
          2],  {e, s}]

... and y/r (see Eq.(6.80)) ...

yor = 
Collect[
        Truncate[Sin[Ω]*cfom + Cos[Ω]*sfom*ci, 
          2],  {e, s}]

... and z/r (see Eq.(6.81)).

zor = 
Collect[
        Truncate[sfom*si, 2],  {e, s}]

There are similar series for x'/r', y'/r' and z'/r'.

xorp = 
        xor  /.  {e -> e', 
            s -> 
              s', ω -> ω', Ω -> Ω', M -> M'};

yorp = 
        yor  /.  {e -> e', 
            s -> 
              s', ω -> ω', Ω -> Ω', M -> M'};

zorp = 
        zor  /.  {e -> e', 
            s -> 
              s', ω -> ω', Ω -> Ω', M -> M'};

This allows us to calculate the series for Cos[ψ] (see Eq.(6.82)).

cpsi = 
      TrigReduce[
        Truncate[xor*xorp + yor*yorp + zor*zorp, 
            2]  /.  {M -> λ - ϖ, ω -> ϖ - Ω, 
            M' -> λ' - ϖ', ω' -> ϖ' - Ω'}]

We can gather up the terms in the eccentricities and inclinations for easier comparison.

Collect[Expand[%],{e,e',s,s'}]

Because θ = ω + Ω + f  we can now find the series for Cos[θ – θ'] (see Eq.(6.83)).

cdtheta=
TrigReduce[Truncate[
(Cos[Ω]*cfom-Sin[Ω]*sfom)*(Cos[Ω']*cfomp-Sin[Ω']*sfomp)+
(Sin[Ω]*cfom+Cos[Ω]*sfom)*(Sin[Ω']*cfomp+Cos[Ω']*sfomp),
2] /. {M->λ-ϖ,ω->ϖ-Ω,M'->λ'-ϖ',ω'->ϖ'-Ω'}];

The quantity Ψ is defined to be the difference between Cos[ψ] and Cos[θ – θ'] (see Eq.(6.84)).

Ψ = TrigReduce[cpsi - cdtheta]

Because this is already second order in the inclinations we can ignore the product (r/a)(r'/a') and write (see Eq.(6.85)):

rapsi2 = Expand[Ψ/2]

Now we need to introduce a series for the equation of centre (f–M).  We start by deriving a series for the eccentric anomaly in terms of the mean anomaly.  (See Eq.(2.80).)

ea=Collect[
    M+Expand[
      2*Normal[
       Series[
        Sum[(1/s)BesselJ[s, s e] Sin[s M],{s,1,5}],
       {e,0,5}]]], e]

Now we need the series for d f  (see Eq.(2.87)).

df=
Collect[TrigReduce[
Normal[Series[Sqrt[1-e^2]*(D[ea, M])^2, {e,0,5}]]],e];

We integrate this term by term to get the series for f in terms of M (see Eq.(2.88)).

f = Truncate[Collect[Integrate[df, M], e], 2]

Express θ in terms of longitudes ...

θ = ω + Ω + 
          f  /.  {M -> λ - ϖ, ω -> ϖ - Ω, 
          M' -> λ' - ϖ', ω' -> ϖ' - Ω'}

... and derive a series for Cos[jθ] ...

cjth = 
      TrigReduce[Normal[Series[Cos[j*θ], {e, 0, 2}]]]

... tidying up the terms (see Eq.(6.88)).

Collect[Expand[%],e]

Similarly, derive a series for Sin[jθ] and tidy it up (see Eq.(6.89)).

sjth = 
      TrigReduce[Normal[Series[Sin[j*θ], {e, 0, 2}]]]

Collect[Expand[%],e]

Derive the series for Cos[jθ'] by substitution.

cjthp = 
      cjth  /.  {e -> e', 
          s -> 
            s', ϖ -> ϖ', Ω -> Ω', λ -> λ'}

Collect[Expand[%],{e,e'}]

Derive the series for Sin[jθ'] by substitution.

sjthp = 
      sjth  /.  {e -> e', 
          s -> 
            s', ϖ -> ϖ', Ω -> Ω', λ -> λ'}

Now we can calculate the series for Cos j[θ-θ'].

cjdth = TrigReduce[Truncate[cjth*cjthp + sjth*sjthp, 2]]

Tidy up the terms (see Eq.(6.90)).

Collect[Expand[%],e]

Now we derive a high-order series for r/a (see Eq.(2.80)).

troa = 
Expand[1+e^2/2-2e*
Series[
Sum[(1/s^2)*
D[BesselJ[s, s e],e]Cos[s*M],{s,1,6}],{e,0,6}]];

Tidy it up ...

Collect[Expand[troa],e]

... and truncate it to second order.

roa = Collect[Expand[Series[Normal[troa],{e,0,2}]], e]

Subtract 1 to get the series for ε (see Eq.(6.91)).

ε = 
      roa - 1  /.  {M -> λ - ϖ}

Derive the series for ε' by substitution.

ε' =  ε  /.  {e -> 
            e', λ -> λ', ϖ -> ϖ'}

Now we are in a position to put all this together to derive an expression for the direct part of the disturbing function to second order.  Here we make use of Eqs.(6.94)–(6.100).  The b[s,j,α,n] are Laplace coefficients.

RD=
(
(1/2)*
(b[1/2,j,α,0]+
ε*α*b[1/2,j,α,1]+
ε'*(-α*b[1/2,j,α,1]-b[1/2,j,α,0])+
(1/2)*ε^2*α^2*b[1/2,j,α,2]+
ε*ε'*(-α^2*b[1/2,j,α,2]-2*α*b[1/2,j,α,1])+
(1/2)*(ε')^2*
(α^2*b[1/2,j,α,2]+4*α*b[1/2,j,α,1]+2*b[1/2,j,α,0])
)+
rapsi2*α*b[3/2,j,α,0])*cjdth;

Truncate and simplify this expression.  N.B. This can take a while!

RD=TrigReduce[Truncate[RD,2]];

Expand the resulting expression for RD (see Eqs.(6.102)–(6.104)).  The remaining work (not shown) consists of identifying and gathering terms in the same argument, using simple substitution if necessary.  This is described in the text.

Expand[RD]

For the indireact parts we need a series for r/a in terms of longitudes (see Eq.(6.108)).

roat=roa /. {M->λ-ϖ}

We also need a series for (((a'/r'))^2) and so we invert the series for r/a square it and substitute (see Eq.(6.109)).

aorp2=TrigReduce[Normal[Series[(1/roat)^2,{e,0,2}]]] /. {e->e',λ->λ',ϖ->ϖ'}

Gather the terms in e' for comparison with Eq.(6.109).

Collect[Expand[%],e']

Put all these together to get the expression for the indirect part with an external perturber (see Eq.(6.110)).

RE=Truncate[-roat*aorp2*cpsi,2]

Re-arrange to gather the zeroth-order arguments.

Collect[RE,Cos[λ-λ']]

Make substitutions to get series for r'/a' and (((a/r))^2)(see Eq.(6.111)).

roatp=roat /. {e->e',λ->λ',ϖ->ϖ'};
aor2=aorp2 /. {e'->e,λ'->λ,ϖ'->ϖ};

Finally we get the series for the indirect part due to an internal perturber (see Eq.(6.111)).

RI=Truncate[-roatp*aor2*cpsi,2]

Re-arrange to gather the zeroth-order arguments.

Collect[RI,Cos[λ-λ']]