Chapter 7, Section 7.8

First we define a function Brouwer that implements the Brouwer & van Woerkom secular solution for the planets (except Pluto) and returns the values of h, k ,e, ϖ, p, q, I, Ω for a given year, t.  Note that Brouwer does not derive a secular theory, it merely implements the published solution of Brouwer & van Woerkom.

Brouwer[t_] :=
Module[
 {g = {5.46326, 7.34474, 17.32832, 18.00233, 4.29591,
       27.77406, 2.71931, 0.63332, -19.18225, 51.25222},
  β = {92.182, 196.881, 335.224, 317.948, 29.550,
          125.120, 131.944, 69.021, 293.979, 220.691},
  f = {-5.21054, -6.57080, -18.74359, -17.63331, 
       0.00000, -25.73355, -2.90266, -0.67752},
  γ = {19.433, 318.057, 255.031, 296.541,
           107.102, 127.367, 315.063, 202.293},
  eij = {{0.17454, 0.00610, 0.00392, 0.00064, -0.00001,
               -0.00001, 0.0, 0.0},
         {-0.02550, 0.02094, 0.01634, 0.00290, -0.00001,
               -0.00001, 0.0, 0.0},
         {0.00154, -0.01255, 0.01043, 0.02972, 0.0,
               -0.00001, 0.0, 0.0},
         {-0.00169, 0.01483, -0.01483, 0.07237, 0.0,
               -0.00001, 0.0, 0.0},
         {0.03571, 0.01909, 0.01834, 0.01861, 0.04482,
               0.03291, -0.03188, 0.00063},
         {0.00096, 0.00006, 0.00283, 0.01502, -0.01535,
               0.04863, -0.00248, -0.00011},
         {0.00043, 0.00042, 0.00044, 0.00059, 0.00147,
               0.00148, 0.03043, -0.00323},
         {0.00001, 0.00001, 0.00001, 0.00002, 0.00005,
               0.00006, 0.00143, 0.00938},
         {-0.00021, -0.00011, -0.00012, -0.00020, -0.00016,
               -0.00049, 0.00020, 0.0},
         {0.00010, 0.00004, 0.00015, 0.00064, -0.00054,
               0.00202, -0.00016, 0.0}},
   iij = {{0.12449, 0.01178, 0.00849, 0.00180, -0.00002, 
               -0.00003, 0.00002, 0.0},
          {-0.03545, 0.01004, 0.00810, 0.00180, -0.00001,
               -0.00002, 0.00001, 0.0},
          {0.00409, -0.02680, 0.02448, -0.03589, 0.0,
               0.0, 0.0, 0.0},
          {0.00116, -0.00685, 0.00453, 0.05025, 0.0,
               -0.00002, 0.0, 0.0},
          {0.02757, 0.02757, 0.02757, 0.02757, 0.02757,
               0.02757, 0.02757, 0.02757},
          {0.00028, 0.00012, 0.00281, 0.00965, -0.00631,
               0.01572, -0.00069, -0.00008},
          {-0.00332, -0.00192, -0.00173, -0.00126, -0.00096,
               -0.00078, 0.01760, -0.00207},
          {-0.00145, -0.00132, -0.00130, -0.00123, -0.00117,
               -0.00113, 0.00110, 0.01175}},
   time
 },
 tg = g; tf = f; tb = β; tga = γ;
 eji = Transpose[eij];
 iji = Transpose[iij];
 time = t-1900.;
 h = eji.(Sin[N[(g*time/3600. + β) Degree]]);
 k = eji.(Cos[N[(g*time/3600. + β) Degree]]);
 e = Table[Sqrt[h[[j]]*h[[j]] + k[[j]]*k[[j]]],{j,8}];
 cp = N[ArcTan[k,h]/Degree];
 p = iji.(Sin[N[(f*time/3600. + γ) Degree]]);
 q = iji.(Cos[N[(f*time/3600. + γ) Degree]]);
 i = Table[N[ArcSin[
       Sqrt[p[[j]]*p[[j]] + q[[j]]*q[[j]]]]/Degree],{j,8}];
 om = N[ArcTan[q,p]/Degree];
 {h,k,e,cp,p,q,i,om}];

Set up a list, time, containing the times we are interested in.  The units here will be in millions of years.

time = Table[ttt,{ttt,-5,5,0.01}];

Create a list, ecdat, of the eccentricity values.

ecdat = Table[Brouwer[tt][[3]],{tt,-5000000,5000000,10000}];

Create a list, idat, of the inclination values.

idat = Table[Brouwer[tt][[7]],{tt,-5000000,5000000,10000}];

Isolate the eccentricities of the planets.

emerc = Transpose[ecdat][[1]];
evenu = Transpose[ecdat][[2]];
eeart = Transpose[ecdat][[3]];
emars = Transpose[ecdat][[4]];
ejupi = Transpose[ecdat][[5]];
esatu = Transpose[ecdat][[6]];
euran = Transpose[ecdat][[7]];
enept = Transpose[ecdat][[8]];

Isolate the inclinations of the planets.

imerc = Transpose[idat][[1]];
ivenu = Transpose[idat][[2]];
ieart = Transpose[idat][[3]];
imars = Transpose[idat][[4]];
ijupi = Transpose[idat][[5]];
isatu = Transpose[idat][[6]];
iuran = Transpose[idat][[7]];
inept = Transpose[idat][[8]];

Plot the eccentricity of Mercury.

pe1 = ListPlot[Transpose[{time,emerc}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Venus.

pe2 = ListPlot[Transpose[{time,evenu}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Earth.

pe3 = ListPlot[Transpose[{time,eeart}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Mars.

pe4 = ListPlot[Transpose[{time,emars}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Jupiter.

pe5 = ListPlot[Transpose[{time,ejupi}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Saturn.

pe6 = ListPlot[Transpose[{time,esatu}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Uranus.

pe7 = ListPlot[Transpose[{time,euran}],
         Frame->True,PlotJoined->True]

Plot the eccentricity of Neptune.

pe8 = ListPlot[Transpose[{time,enept}],
         Frame->True,PlotJoined->True]

Plot the inclination of Mercury.

pi1 = ListPlot[Transpose[{time,imerc}],
         Frame->True,PlotJoined->True]

Plot the inclination of Venus.

pi2 = ListPlot[Transpose[{time,ivenu}],
         Frame->True,PlotJoined->True]

Plot the inclination of Earth.

pi3 = ListPlot[Transpose[{time,ieart}],
         Frame->True,PlotJoined->True]

Plot the inclination of Mars.

pi4 = ListPlot[Transpose[{time,imars}],
         Frame->True,PlotJoined->True]

Plot the inclination of Jupiter.

pi5 = ListPlot[Transpose[{time,ijupi}],
         Frame->True,PlotJoined->True]

Plot the inclination of Saturn.

pi6 = ListPlot[Transpose[{time,isatu}],
         Frame->True,PlotJoined->True]

Plot the inclination of Uranus.

pi7 = ListPlot[Transpose[{time,iuran}],
         Frame->True,PlotJoined->True]

Plot the inclination of Neptune.

pi8 = ListPlot[Transpose[{time,inept}],
         Frame->True,PlotJoined->True]

Now we can consider the Jupiter-Saturn system under the Brouwer & van Woerkom solution.

Set up a list of the times we are interested in.  Here we choose the same interval as we used in Fig.7.1.

jstime = Table[ttt,{ttt,-0.1,0.1,0.001}];

Calculate the eccentricities ...

ecjsdat = Table[Brouwer[tt][[3]],{tt,-100000,100000,1000}];

... and inclinations at these times.

ijsdat = Table[Brouwer[tt][[7]],{tt,-100000,100000,1000}];

Isolate Jupiter's eccentricity.

ejupi = Transpose[ecjsdat][[5]];

Isolate Saturn's eccentricity.

esatu = Transpose[ecjsdat][[6]];

Isolate Jupiter's inclination.

ijupi = Transpose[ijsdat][[5]];

Isolate Saturn's inclination.

isatu = Transpose[ijsdat][[6]];

Derive the plot for Jupiter's eccentricity but don't display it.

pej=ListPlot[Transpose[{jstime,ejupi}],
Frame->True,PlotJoined->True,DisplayFunction->Identity]

Derive the plot for Saturn's eccentricity but don't display it.

pes=ListPlot[Transpose[{jstime,esatu}],
Frame->True,PlotJoined->True,DisplayFunction->Identity]

Display the eccentricity plot.

Show[pej,pes,DisplayFunction->$DisplayFunction]

Derive the plot for Jupiter's inclination but don't display it.

pij=ListPlot[Transpose[{jstime,ijupi}],
Frame->True,PlotJoined->True,DisplayFunction->Identity]

Derive the plot for Saturn's inclination but don't display it.

pis=ListPlot[Transpose[{jstime,isatu}],
Frame->True,PlotJoined->True,DisplayFunction->Identity]

Display the inclination plot.

Show[pij,pis,DisplayFunction->$DisplayFunction]