Chapter 7, Section 7.10, Figure 7.12

We start by defining laplace, a procedure to calculate a Laplace coefficient by evaluating an integral (see Eq.(6.67)).

laplace[ s_, jj_, aa_ ] := 
    If[aa != 0,
       NIntegrate[ N[ (2/Pi)*(Cos[jj*x])/((1-2aa*Cos[x]+aa^2)^s) ], 
               {x, 0, Pi} ],
       0];

Now we define a function BrouwerbigAA, a slightly shorter version of the function Brouwer from Sect.7.8.nb, that returns the free precession rate for the value atest of the test particle's semi-major axis.

BrouwerbigAA[atest_, year_] :=
Module[
 {m = 1/{6025607.0, 408662.0, 328912.0, 3099774.0, 1047.91,
         3500.05, 22970.0, 19360.0},
  a = {0.38710, 0.72333, 1.00002, 1.52364, 5.20248,
       9.56329, 19.29373, 30.27435},
  g = {5.46326, 7.34474, 17.32832, 18.00233, 4.29591,
       27.77406, 2.71931, 0.63332, -19.18225, 51.25222},
  f = {-5.21054, -6.57080, -18.74359, -17.63331, 
       0.00000, -25.73355, -2.90266, -0.67752},
   ntest=360./atest^(3/2)
},
tg = g; tf = f; tb = beta; tga = gamma;
gr = N[g Degree];
fr = N[f Degree];
time = year-1900.;
alphaj    = (atest/a /. {x_ :> 1/x /; x>1});
alphabarj = (atest/a /. {x_ :> 1   /; x>1});
b321 = Map[laplace[3/2, 1, #]&, alphaj];
bigAA =  ntest*((1/4)*(Plus @@ (m*alphaj*alphabarj*b321)));
asbigAA = bigAA*60*60;
N[{atest, asbigAA}]
];

Calculate the value of the free precession rate between 0.1 and 40 AU for epoch 1997.

precess = 
Table[BrouwerbigAA[avalue, 1997], {avalue, 0.1, 40, 0.1}];

Plot the resulting values in the range 4 to 40 AU.

pout =
ListPlot[precess, PlotJoined->True, 
Frame->True, PlotRange->{{4,40},{0,60}}]

Set the secular eigenfrequencies of the Brouwer and van Woerkom system.

gg = {5.46326, 7.34474, 17.32832, 18.00233, 4.29591,
       27.77406, 2.71931, 0.63332, -19.18225, 51.25222};
ff = {-5.21054, -6.57080, -18.74359, -17.63331, 
       0.00000, -25.73355, -2.90266, -0.67752};

Create a graphics object of lines denoting the secular eigenfrequencies.

gout = Graphics[Map[Line[{{4,#},{40,#}}]&, gg]];
fout = Graphics[
Join[{Dashing[{0.01, 0.01}]},
     Map[Line[{{4,#},{40,#}}]&, -ff]]];

Display the plot with the eigenfrequencies.

Show[pout, gout, fout]

Now consider the inner part of the solar system.  Calculate the free precession rate between 0.1 and 4 AU.

precessinner = 
Table[BrouwerbigAA[avalue, 1997], {avalue, 0.1, 4, 0.01}];

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

pin = 
ListPlot[precessinner, PlotJoined->True, 
Frame->True, PlotRange->{{0,4},{0,60}}]

Create a graphics object of lines denoting the secular eigenfrequencies.

gin = Graphics[Map[Line[{{0,#},{4,#}}]&, gg]];
fin = Graphics[
Join[{Dashing[{0.01, 0.01}]},
     Map[Line[{{0,#},{4,#}}]&, -ff]]];

Display the plot with the eigenfrequencies.

Show[pin, gin, fin]