Chapter 7, Section 7.10, Figure 7.16

We start by defining laplace, a function 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 BrouwerProper, a slightly longer version of the function Brouwer from Sect.7.8.nb, that returns the forced orbital elements of a test particle with semi-major axis atest at year year using the Brouwer & van Woerkom secular theory for the planets.

BrouwerProper[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},
  β = {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}},
   ntest=360./atest^(3/2),time
},
tg = g; tf = f; tb = β; tga = γ;
gr = N[g Degree];
fr = N[f Degree];
eji = Transpose[eij];
iji = Transpose[iij];
time = year-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];
alphaj    = (atest/a /. {x_ :> 1/x /; x>1});
alphabarj = (atest/a /. {x_ :> 1   /; x>1});
b321 = Map[laplace[3/2, 1, #]&, alphaj];
b322 = Map[laplace[3/2, 2, #]&, alphaj];
bigAA =  ntest*((1/4)*(Plus @@ (m*alphaj*alphabarj*b321)));
bigBB = -ntest*((1/4)*(Plus @@ (m*alphaj*alphabarj*b321)));
bigAj = -(1/4)*m*ntest*alphaj*alphabarj*b322;
bigBj =  (1/4)*m*ntest*alphaj*alphabarj*b321;
nui = bigAj . eji;
mui = bigBj . iji;
h0 = -(Plus @@ ((nui/(bigAA-g/3600.))*
       Sin[N[(g*time/3600. + β ) Degree]]));
k0 = -(Plus @@ ((nui/(bigAA-g/3600.))*
       Cos[N[(g*time/3600. + β ) Degree]]));
eforced = Sqrt[h0^2+k0^2];
periforced = N[ArcTan[k0,h0]/Degree];
p0 = -(Plus @@ ((mui/(bigBB-f/3600.))*
       Sin[N[(f*time/3600. + γ) Degree]]));
q0 = -(Plus @@ ((mui/(bigBB-f/3600.))*
       Cos[N[(f*time/3600. + γ) Degree]]));
iforced = N[ArcSin[Sqrt[p0^2+q0^2]]/Degree];
nodeforced = N[ArcTan[q0,p0]/Degree];
N[{atest, eforced, periforced, iforced, nodeforced}]
];

Calculate the forced elements between 2 and 4 AU for the epoch of 1997.

results = 
N[Table[BrouwerProper[aaa, 1997],{aaa, 2, 4, 0.01}]];

Transpose the results ...

tres = Transpose[results];

... and prepare lists of forced eccentricity, forced longitude of pericentre, forced inclination and forced longitude of ascending node as a function of semi-major axis.

ae = Transpose[{tres[[1]], tres[[2]]}];
acp = Transpose[{tres[[1]], tres[[3]]}];
ai = Transpose[{tres[[1]], tres[[4]]}];
aom = Transpose[{tres[[1]], tres[[5]]}];

Select those asteroids in the first 64 and then the remaining parts of ae with forced eccentricities less than 0.09.  The purpose of breaking things up is to take account of the singularity caused by the secular resonance near 2.6 AU.

ae1 = Select[Take[ae,64], (#[[2]]<0.09)&];

ae2 = Select[Drop[ae,64], (#[[2]]<0.09)&];

Obtain the corresponding set of longitudes of pericentre.

acp1 = Take[acp,64];

acp2 = Drop[acp,64];

Plot the forced eccentricity as a function of semi-major axis for the first 64 ...

pe1 = 
ListPlot[ae1, PlotRange->{{2,4},{0,0.08}},
Frame->True, PlotJoined->True, AxesOrigin->{2,0}]

... and the remainder

pe2 = 
ListPlot[ae2, PlotRange->{{2,4},{0,0.08}},
Frame->True, PlotJoined->True, AxesOrigin->{2,0}]

Show both plots together.

pe = Show[pe1,pe2]

Now do the same for the forced longitudes of pericentre.

pcp1 = 
ListPlot[acp1, PlotRange->{{2,4},{-40,60}},
Frame->True, PlotJoined->True, AxesOrigin->{2,-40}]

pcp2 = 
ListPlot[acp2, PlotRange->{{2,4},{-40,60}},
Frame->True, PlotJoined->True, AxesOrigin->{2,-40}]

Show both plots together.

pcp = Show[pcp1, pcp2]

The secular resonance does not affect the forced inclination and node and so these plots are straight-forward.

pi = 
ListPlot[ai, PlotRange->{{2,4},{0,3}}, 
PlotJoined->True, Frame->True]

pom = 
ListPlot[aom, PlotRange->{{2,4},{-40,110}}, 
PlotJoined->True, Frame->True, AxesOrigin->{2,-40}]