Chapter 7, Section 7.3

We start by defining a function secular that derives a secular perturbation theory for any number of masses with given orbital elements.

(* 
               ==== Secular perturbation solver ==== 
                       v1.0 February 1993 
                         Carl D. Murray             
*)


secular::usage = "secular[nbody, m, a, ecc, inc, peri, node, n,
	rp, j2, j4] determines the secular perturbation solution for
	the gravitational interaction of nbody bodies with mass ratios
	m, semi-major axes, a, eccentricities, ecc, inclinations, inc,
	longitudes of pericentre, peri, longitudes of ascending node,
	node, mean motions, n (in degrees/year) which are orbiting a
	central body of radius, rp with J2 and J4 values, j2 and j4.
	The function defines ecalc and icalc which can be used to plot
	the time variation of the eccentricity and inclination of
	any of the bodies.";

secular[nbody_, m_, a_, ecc_, inc_, peri_, node_, n_, 
        rp_, j2_, j4_] :=
(* 
   Based on the theory given in Dermott & Nicholson, Nature,
   319, 115-120.
   
   nbody : number of bodies
   m     : masses of bodies
   a     : semi-major axes of bodies
   ecc   : eccentricities of bodies
   inc   : inclinations of bodies
   peri  : longitudes of pericentre of bodies
   node  : longitudes of ascending nodes of bodies
   n     : mean motion (in degrees/year) of bodies
   rp    : radius of central object
   j2    : J2 of central object
   j4    : J4 of central object
*)  
Module[{rpa2 = N[(rp/a)^2], rpa4 = N[(rp/a)^4], j22=j2*j2,
        mred = 1+m},
Print["Mass ratios :"];
Print[m];
Print[" "];
Print["Semi-major axes :"];
Print[a];
Print[" "];
Print["Mean motions (degrees/year) :"];
Print[n];
Print[" "];
Print["Initial values :"];
Print["==============  "];
Print[" "];
Print["Eccentricities :"];
Print[ecc];
Print[" "];
Print["Longitudes of pericentre :"];
Print[peri];
Print[" "];
Print["Inclinations :"];
Print[inc];
Print[" "];
Print["Longitudes of ascending node :"];
Print[node];
Print[" "];
(* 
   Set the initial values of h, k, p, q
*)
hinit = ecc*Sin[N[peri]];
kinit = ecc*Cos[N[peri]];
pinit = N[inc*Sin[node]];
qinit = N[inc*Cos[node]];
(*
   Define the alpha, alphabar functions from the D&N paper
*)
alpha[ j_, k_ ] := a[[k]]/a[[j]] /; k <  j;
alpha[ j_, k_ ] := 0             /; k == j;
alpha[ j_, k_ ] := a[[j]]/a[[k]] /; k >  j;
alphabar[ j_, k_ ] := 1             /; k <  j;
alphabar[ j_, k_ ] := 0             /; k == j;
alphabar[ j_, k_ ] := alpha[ j, k ] /; k >  j;
(*
  Define the Laplace coefficient function using a numerical
  integration
*)
laplace[ s_, jj_, aa_ ] := 
  If[aa != 0,
     NIntegrate[ N[ (1/Pi)*(Cos[jj*x])/((1-2aa*Cos[x]+aa^2)^s) ],
                {x, 0, 2Pi} ],
     0];
(*
  Define the elements of the matrices A and B
*)
A[ j_, k_ ] := -(1/4)*m[[k]]*(n[[j]]/mred[[j]])*
                 alpha[j,k]*alphabar[j,k]*
                 laplace[3/2, 2, alpha[j,k]] /; j != k;
A[ j_, k_ ] := n[[j]]*
               ((3/2)*j2*rpa2[[j]] - (9/8)*j22*rpa4[[j]] -
                (15/4)*j4*rpa4[[j]] +
                ((1/4)/mred[[j]])*
                         Sum[If[l != j,
                            m[[l]]*
                            alpha[j,l]*alphabar[j,l]*
                            laplace[3/2, 1, alpha[j,l]],
                            0],
                          {l,nbody}]) /; j == k;
B[ j_, k_ ] := (1/4)*m[[k]]*(n[[j]]/mred[[j]])*
                 alpha[j,k]*alphabar[j,k]*
                 laplace[3/2, 1, alpha[j,k]] /; j != k;
B[ j_, k_ ] := -n[[j]]*
               ((3/2)*j2*rpa2[[j]] - (27/8)*j22*rpa4[[j]] -
                (15/4)*j4*rpa4[[j]] +
                ((1/4)/mred[[j]])*
                         Sum[If[l != j,
                            m[[l]]*
                            alpha[j,l]*alphabar[j,l]*
                            laplace[3/2,1,alpha[j,l]],
                            0],
                          {l,nbody}]) /; j == k;
(*
  Construct the two matrices
*)
bigA = Array[ A, {nbody,nbody} ];
bigB = Array[ B, {nbody,nbody} ];
Print["Matrix A (deg/yr):"];
Print[MatrixForm[bigA]];
Print[" "];
Print["Matrix B (deg/yr):"];
Print[MatrixForm[bigB]];
Print[" "];
(*
  Solve the resulting eigensystems
*)
esA = Eigensystem[ bigA ];
esB = Eigensystem[ bigB ];
(*
  Extract the two sets of eigenfrequencies ...
*)
g = esA[[1]];
f = esB[[1]];
Print["g eigenfrequencies (deg/yr) :"];
Print[g];
Print[" "];
Print["g eigenvectors :"];
Print[MatrixForm[esA[[2]]]];
Print[" "];
Print["f eigenfrequencies (deg/yr) :"];
Print[f];
Print[" "];
Print["f eigenvectors :"];
Print[MatrixForm[esB[[2]]]];
Print[" "];
Print["Initial h values :"];
Print[hinit];
Print[" "];
Print["Initial k values :"];
Print[kinit];
Print[" "];
Print["Initial p values :"];
Print[pinit];
Print[" "];
Print["Initial q values :"];
Print[qinit];
Print[" "];
(*
  ... and the corresponding (unscaled) eigenvectors
*)
ebarji = Transpose[esA[[2]]];
ibarji = Transpose[esB[[2]]];
(*
  Determine the values of Sin[beta] and Sin[gamma] and hence
  the scaling of the eigenvectors
*)
fsbeta = LinearSolve[ebarji, hinit];
fcbeta = LinearSolve[ebarji, kinit];
beta = ArcTan[fcbeta, fsbeta];
Print["Beta (deg) :"];
Print[N[beta/Degree]];
Print[" "];
sbeta = Sin[beta];
fsgamma = LinearSolve[ibarji, pinit];
fcgamma = LinearSolve[ibarji, qinit];
gamma = ArcTan[fcgamma, fsgamma];
Print["Gamma (deg) :"];
Print[N[gamma/Degree]];
Print[" "];
sgamma = Sin[gamma];
emag = fsbeta/sbeta;
imag = fsgamma/sgamma;
eji = Map[Times[#,emag]&, ebarji];
iji = Map[Times[#,imag]&, ibarji];
gr = N[g Degree];
fr = N[f Degree];
Print["Scaled h-k eigenvectors :"];
Print[MatrixForm[eji]];
Print[" "];
Print["Scaled p-q eigenvectors :"];
Print[MatrixForm[iji]];
Print[" "];
(*
  Use these values to determine e and I as functions of time
*)
ejcalc[j_, t_] := 
  Sqrt[(Plus @@ (eji[[j]]*Sin[gr*t+beta]))^2 + 
       (Plus @@ (eji[[j]]*Cos[gr*t+beta]))^2];
ijcalc[j_, t_] := 
  N[
   Sqrt[(Plus @@ (iji[[j]]*Sin[fr*t+gamma]))^2 + 
        (Plus @@ (iji[[j]]*Cos[fr*t+gamma]))^2]]
];

Now apply the function secular to the Jupiter-Saturn system.  Note that we are using mean motions derived from Kepler's third law and we are setting J2 and J4 to be zero (i.e. there is no oblateness of the central body to consider).

secular[2, 
        {1./(1047.355), 1./(3498.500)},
        {5.202545297, 9.554840607},
        {0.047462192, 0.057548133},
        {1.30667, 2.48795} Degree,
        {13.983865305, 88.719424868} Degree,
        {100.0381, 113.1334} Degree,
        {360./5.202545297^(3/2), 360./9.554840607^(3/2)},
        0,
        0,
        0];

The function secular has also defined a function ejcalc that returns the eccentricity of any of the bodies at a specific time. We can use this to derive Eq.(7.52).

Chop[FullSimplify[ejcalc[1,t]]]

Chop[FullSimplify[ejcalc[2,t]]]

Convert some of the constants for a proper comparison.

N[0.0000896338/Degree]

N[1.63251/Degree]

The function secular has also defined a function ijcalc that returns the inclination of any of the bodies at a specific time. We can use this to derive Eq.(7.53).

Chop[FullSimplify[ijcalc[1,t]]]

Chop[FullSimplify[ijcalc[2,t]]]

Convert some of the constants for a proper comparison.

N[0.368006/Degree]

N[0.000123264/Degree]

Now plot the evolution of the eccentricities of Jupiter and Saturn as a function of time (see Fig.7.1a).

Plot[{ejcalc[1,t],ejcalc[2,t]},{t,-100000,100000},
      PlotRange->{{-100000,100000},{0,0.1}},Frame->True]

Do the same for the inclinations of Jupiter and Saturn (see Fig.7.1b).

Plot[{ijcalc[1,t]/Degree,ijcalc[2,t]/Degree},{t,-100000,100000},
      PlotRange->{{-100000,100000},{0,3}},Frame->True]