Chapter 7, Section 7.5

We start by defining a function secular that derives a secular perturbation theory for any number of masses with given orbital elements.  In this version we have removed all the printing commands.  Sect.7.3.nb contains the version that prints out the results.

(* 
               ==== 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},
(* 
   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} ];
(*
  Solve the resulting eigensystems
*)
esA = Eigensystem[ bigA ];
esB = Eigensystem[ bigB ];
(*
  Extract the two sets of eigenfrequencies ...
*)
g = esA[[1]];
f = esB[[1]];
(*
  ... 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];
sbeta = Sin[beta];
fsgamma = LinearSolve[ibarji, pinit];
fcgamma = LinearSolve[ibarji, qinit];
gamma = ArcTan[fcgamma, fsgamma];
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];
(*
  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]/Degree]
];

Apply secular to derive a seculr theory for the Jupiter-Saturn system.

secular[2, 
        {1./(1047.355), 1./(3498.500)},
        {5.202545, 9.554841},
        {0.0474622, 0.0575481},
        {1.30667, 2.48795} Degree,
        {13.983865, 88.719425} Degree,
        {100.0381, 113.1334} Degree,
        {30.3374, 12.1890},
        0,
        0,
        0
        ]

Now define a function proper that will caculate proper elements for a test particle in this system with a given set of orbital elements.

proper[nbody_, m_, a_, ecc_, inc_, peri_, node_, n_, rp_, j2_, j4_,
       atest_, etest_, itest_, ptest_, nodtest_, dummy_, time_] :=
Module[{ntest=360./atest^(3/2)},
(*
secular[nbody, m, a, ecc, inc, peri, node, n, rp, j2, j4];
*)
alphaj =    (atest/a /. {x_ :> 1/x /; x>1});
alphabarj = (atest/a /. {x_ :> 1   /; x>1});
bigNj = (1/8)*alphaj*alphabarj*Map[laplace[3/2, 1, #]&, alphaj];
bigPj = (1/8)*alphaj*alphabarj*Map[laplace[3/2, 2, #]&, alphaj];
bigAA = (2*ntest) * (Plus @@ (m*bigNj));
nuj   = (2*ntest) * ( (m*bigPj) . (eji) );
muj   = (2*ntest) * ( (m*bigNj) . (iji) );
h0 = Plus @@ ( (nuj/(bigAA-g)) Sin[gr*time+beta] );
k0 = Plus @@ ( (nuj/(bigAA-g)) Cos[gr*time+beta] );
eforced = Sqrt[h0^2+k0^2];
pforced = N[ArcTan[k0, h0]/Degree];
p0 = Plus @@ ( (muj/(bigAA+f)) Sin[fr*time+gamma] );
q0 = Plus @@ ( (muj/(bigAA+f)) Cos[fr*time+gamma] );
iforced = N[Sqrt[p0^2+q0^2]/Degree];
nodforced = N[ArcTan[q0, p0]/Degree];
h = etest*Sin[N[ptest]];
k = etest*Cos[N[ptest]];
p = N[itest*Sin[nodtest]];
q = N[itest*Cos[nodtest]];
nusin = h - h0;
nucos = k - k0;
efree = Sqrt[nusin^2+nucos^2];
musin = p - p0;
mucos = q - q0;
ifree = N[ArcSin[Sqrt[musin^2+mucos^2]]/Degree];
pbeta  = N[ArcTan[nucos, nusin]/Degree]-bigAA*time;
pgamma = N[ArcTan[mucos, musin]/Degree]+bigAA*time;
pbetar  = N[pbeta  Degree];
pgammar = N[pgamma Degree];
{atest,bigAA,eforced,pforced,iforced,nodforced}
];

Create a list, data, of the values of the semi-major axis (atest), the free precession rate (bigAA), the forced eccentricity (eforced), pericentre (pforced), inclination (iforced) and node (nodforced) for semi-major axes up to 30 AU.

data = Table[
       proper[2, 
        {1./(1047.355), 1./(3498.500)},
        {5.202545, 9.554841},
        {0.0474622, 0.0575481},
        {1.30667, 2.48795} Degree,
        {13.983865, 88.719425} Degree,
        {100.0381, 113.1334} Degree,
        {30.3374, 12.1890},
        0,
        0,
        0,
        ai, 0.15, 10 Degree, 42 Degree, -10 Degree, 149.07, 0.],
{ai,0.0001,30.,0.05}];

Take the transpose of the resulting file ...

tdata = Transpose[data];

... and extract the values of the semi-major axis, free precession rate, and the forced values of the eccentricity, longiutde of pericentre, inclination and longiutde of ascending node.

adat = tdata[[1]];
bigAdat = tdata[[2]];
efdat = tdata[[3]];
pfdat = tdata[[4]];
ifdat = tdata[[5]];
nfdat = tdata[[6]];

Plot the free precession rate as a function of semi-major axis.

bigAplot = ListPlot[Transpose[{adat,bigAdat}],
             Frame->True,PlotRange->{{0,30},{0,0.05}},
             PlotJoined->True];

Create a graphics object to denote the eigenfrequencies of the system.

gflines = Graphics[
           {Line[{{0.0,Abs[g[[1]]]},{30.0,Abs[g[[1]]]}}],
            Line[{{0.0,Abs[g[[2]]]},{30.0,Abs[g[[2]]]}}],
            Line[{{0.0,Abs[f[[1]]]},{30.0,Abs[f[[1]]]}}],
            Line[{{0.0,Abs[f[[1]]]},{30.0,Abs[f[[1]]]}}]}];

Combine the free precession arte plots with the eigenfrequencies (see Fig.7.4).

Show[bigAplot,gflines];

Plot the forced eccentricity as a function of semi-major axis (see Fig.7.5a).

efplot = ListPlot[Transpose[{adat,efdat}],
             Frame->True,PlotRange->{{0,30},{0,0.1}},
             PlotJoined->True];

Plot the forced longitude of pericentre as a function of semi-major axis (see Fig.7.5b).

pfplot = ListPlot[Transpose[{adat,Mod[pfdat,360]}],
             Frame->True,PlotRange->{{0,30},{0,360}},
             PlotJoined->True];

Plot the forced inclination as a function of semi-major axis (see Fig.7.6a).

ifplot = ListPlot[Transpose[{adat,ifdat}],
             Frame->True,PlotRange->{{0,30},{0,3}},
             PlotJoined->True];

Plot the forced longitude of ascending node as a function of semi-major axis (see Fig.7.6b).

nfplot = ListPlot[Transpose[{adat,Mod[nfdat,360]}],
             Frame->True,PlotRange->{{0,30},{0,360}},
             PlotJoined->True];