Chapter 6, Section 6.11

We start by deriving explicit expressions for the Legendre polynomials of degree 2 in terms of Sin[α] (see Eq.(6.252)) ...

p2 = LegendreP[2,Sin[α]]

... and degree 4 (see Eq.(6.253)).

p4 = LegendreP[4,Sin[α]]

Use the substitution given in Eq.(6.254) and expand the resulting expressions.

ps2 = Expand[p2 /. {Sin[α]->si*(sf*cw+cf*sw)}]

ps4 = Expand[p4 /. {Sin[α]->si*(sf*cw+cf*sw)}]

Use Eq.(6.244) to derive an expression for 𝒢(Subscript[m, p])

gmp = 
Normal[Series[
n^2*a^3/(1+(3/2)*J2*(Rpoa)^2-
 (15/8)*J4*(Rpoa)^4),{Rpoa,0,4}]] /. {Rpoa->Rp/a}

Now use these expressions to derive an expression for ℛ as defined in Eq.(6.251).

bigR = 
-GM*(J2*Rp^2/r^3*ps2 +J4*Rp^4/r^5*ps4) /. 
{r->a(1-e^2)/(1+e*cf), GM->gmp};

Substitute the series for Sin[f] and Cos[f] (to second order in the eccentricity) as given by Eqs.(2.84) and (2.85).

bigR = bigR /. 
{sf->Sin[M]+e*Sin[2M]+e^2((9/8)*Sin[3M]-(7/8)*Sin[M]),
 cf->Cos[M]+e*(Cos[2M]-1)+(9/8)*e^2*(Cos[3M]-Cos[M])};

Expand the result as a series in e and Sin[I] to second order.

tbigR = 
Normal[Series[(bigR /. {e->t*e, si->t*si}), {t,0,2}]] /. {t->1};

Now expand this as a series in (Subscript[R, p])/a.

tbigR=TrigReduce[Normal[Series[tbigR,{Rp,0,4}]]]

Now that we have a series in the mean anomaly we can carry out a time average of this expression.  Note that the result has terms in (e^2)and (s^2)as well as terms which are independent of e and s and are effectively constant

avbigR = 
Expand[Integrate[tbigR/2/Pi, {M,0,2Pi}] /. 
{cw^2->1-sw^2}]

Isolate the terms in eccentricity  (see the (e^2)terms in Eq.(6.255)) ...

eterm = Coefficient[avbigR, e, 2]

... and inclination  (see the (s^2)terms in Eq.(6.255)).

iterm = Coefficient[avbigR, si, 2]

The remaining terms are the "constant" part.

Expand[final - e^2*eterm - si^2*iterm]

Our expressions for eterm and iterm are both in terms of n, the actual mean motion (including the effects of oblateness) rather than n0, the nominal, keplerian mean motion.  So we use the relationaship given in Eq.(6.244) to express eterm as a series in terms of n0.

Normal[Series[
eterm /. {n^2->n0^2*(1+(3/2)*J2*(Rp/a)^2-(15/8)*J4*(Rp/a)^4)},
{Rp,0,4}]]

Now do the same for iterm.

Normal[Series[
iterm /. {n^2->n0^2*(1+(3/2)*J2*(Rp/a)^2-(15/8)*J4*(Rp/a)^4)},
{Rp,0,4}]]

Input the equations for (n^2), (κ^2) and (μ^2)given in Eqs.(6.244) to (6.246).

n2 = 
  (GM/a^3)*(1+(3/2)J2*(Rp/a)^2-(15/8)*J4*(Rp/a)^4);
ka2 =
  (GM/a^3)*(1-(3/2)J2*(Rp/a)^2+(45/8)*J4*(Rp/a)^4);
mu2 =
  (GM/a^3)*(1+(9/2)J2*(Rp/a)^2-(75/8)*J4*(Rp/a)^4);

Now use Eq.(6.247) to derive the equation for the variation of the longitude of pericentre given in Eq.(6.256).

ϖdot = 
Expand[Normal[Series[
Sqrt[n2]-Sqrt[ka2], {Rp,0,4}]]] /. GM->gmp;

Ensure that (Sqrt[(n^2)])is taken to be n  to get Eq.(6.256).

Normal[Series[ϖdot, {Rp,0,4}]] /. Sqrt[n^2]->n

Now use Eq.(6.248) to derive the equation for the variation of the longitude of ascending node given in Eq.(6.257).

Ωdot = 
Expand[Normal[Series[
Sqrt[n2]-Sqrt[mu2], {Rp,0,4}]]] /. GM->gmp;

Ensure that (Sqrt[(n^2)])is taken to be n  to get Eq.(6.257).

Normal[Series[Ωdot, {Rp,0,4}]] /. Sqrt[n^2]->n