Chapter 8, Section 8.7, Figure 8.7

Input the data on first-order resonances from Table 8.5.  For each resonance we input the coefficient of λ', the coefficient of λ, the value of α = a/a', the value of α(Subscript[f, d])(α), and the maximum value of the eccentricity in the plot (0.3 in this case).

first = 
{{2, -1, 0.629961, -0.749964, 0.3}, 
 {3, -2, 0.763143, -1.54553,  0.3},
 {4, -3, 0.825482, -2.34472,  0.3}};

TableForm[first]

Now do the same for the second-order resonances in which we are interested.

second = 
{{3, -1, 0.480740, 0.287852, 0.3},
 {5, -3, 0.711379, 2.32892,  0.3}};

TableForm[second]

Define a function, FirstPlot, that will do the plots of maximum libration width (based on the pendulum model) in the (a, e) plane.

FirstPlot[mass_, {j1_, j2_, α_, afa_, emax_}] :=
Module[{crn = Abs[mass*afa], evals, p1, p2, a=5.20336*α},
evals = Table[ecc, {ecc, 0.001, emax, 0.001}];
plus = 
N[Map[(a + a*(-(2/9)*crn/j2/# +
Sqrt[(16/3)*crn*#]*Sqrt[1+crn/27/j2/j2/#^3]))&, evals]];
minus = 
N[Map[(a + a*(-(2/9)*crn/j2/# -
Sqrt[(16/3)*crn*#]*Sqrt[1+crn/27/j2/j2/#^3]))&, evals]];
p1 = ListPlot[Transpose[{plus,  evals}], 
      PlotJoined->True, DisplayFunction->Identity];
p2 = ListPlot[Transpose[{minus, evals}], 
      PlotJoined->True, DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]
]

Do the plot for the 2:1 resonance.  The first parameter is the mass ratio of the perturber (Jupiter).

p21 = FirstPlot[0.000954, first[[1]]]

Similarly do the plot for the 3:2 resonance ...

p32 = FirstPlot[0.000954, first[[2]]]

... and the 4:3 resonance.

p43 = FirstPlot[0.000954, first[[3]]]

Now define a function, SecondPlot, to do the same for a second-order resonance.

SecondPlot[mass_, {j1_, j2_, al_, afa_, emax_}] :=
Module[{crn = Abs[mass*afa], evals, p1, p2, a=5.20336*al},
evals = Table[ecc, {ecc, 0.001, emax, 0.001}];
plus  = N[Map[(a + a*Sqrt[(16/3)*crn*#^2])&, evals]];
minus = N[Map[(a - a*Sqrt[(16/3)*crn*#^2])&, evals]];
p1 = ListPlot[Transpose[{plus,  evals}], 
      PlotJoined->True, DisplayFunction->Identity];
p2 = ListPlot[Transpose[{minus, evals}], 
      PlotJoined->True, DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]
]

Plot the maximum libration width of the 3:1 resonance ...

p31 = SecondPlot[0.000954, second[[1]]]

... and the 5:3 resonance.

p53 = SecondPlot[0.000954, second[[2]]]

Put all the plots together in one diagram.

Show[p31,p21,p53,p32,p43,PlotRange->{{2.4,4.7},{0,0.3}},
Frame->True, AxesOrigin->{0.45,0}]