Chapter 8, Section 8.2, Figure 8.4a,c

We start by defining the radius as a function of semi-major axis, eccentricity and true anomaly.

r[a_,e_,t_]:=a(1-e^2)/(1+e Cos[t]);

Define a graphics object to represent the Sun.

sun=Graphics[{PointSize[0.05],Point[{0,0}]}];

Now we define a function ksolve that solves Kepler's equation to find the eccentrci anomaly for a given value of the mean anomaly and eccentricity.

ksolve[ma_, e_] := FindRoot[ea - e Sin[ea] == ma, {ea,ma}];

The function position calculates the x and y coordinates of an object for a particular value of the mean anomaly.

position[ma_,a_,e_]:=
(eca=ea /. ksolve[ma,e];
x=N[a(Cos[eca]-e)]; y=N[a Sqrt[1-e^2] Sin[eca]];
{x,y})

Define a function PlotPath to show points along the orbital path in the rotating frame of the resonant orbit defined by the list of two numbers in pl.  The points are plotted at time intervals of dt.

PlotPath[pl_List, a_, e_, dt_] :=
Module[{time, p1, p2, p = Apply[Rational,pl],
        pt = If[pl[[1]] > pl[[2]], pl[[1]], pl[[2]]],
        sun=Graphics[{PointSize[0.05],Point[{0,0}]}],
        box = 
        Graphics[Line[{{-1.5,-1.5},{1.5,-1.5},
                       {1.5,1.5},{-1.5,1.5},{-1.5,-1.5}}]]},
  time = Table[N[t], {t,0,pt*2*Pi,dt}];
  posvec = Map[position[#,a,e]&,time];
  f = Map[N[ArcTan[#[[1]],#[[2]]]]&,posvec];
  rad = Map[N[Sqrt[#[[1]]^2+#[[2]]^2]]&,posvec];
  points = Transpose[{rad*Cos[p*time-f],rad*Sin[p*time-f]}];
  p1 = 
  ListPlot[points,
           AspectRatio->Automatic, Axes->None,
           Prolog->AbsolutePointSize[5], 
           DisplayFunction->Identity,
           PlotRange->{{-1.5,1.5},{-1.5,1.5}}];
  Show[box,sun,p1,DisplayFunction->$DisplayFunction,
  AspectRatio->Automatic,PlotRange->All]]

Plot points along the paths of orbits at the 2:1internal resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{1,2}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{1,2}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{1,2}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{1,2}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 3:2 internal resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{2,3}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{2,3}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{2,3}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{2,3}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 4:3 internal resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{3,4}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{3,4}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{3,4}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{3,4}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 5:4 internal resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{4,5}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{4,5}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{4,5}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{4,5}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 6:5 internal resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{5,6}, 1.0, 0.1, Pi/15];
p2 = PlotPath[{5,6}, 1.0, 0.2, Pi/15];
p3 = PlotPath[{5,6}, 1.0, 0.3, Pi/15];
p4 = PlotPath[{5,6}, 1.0, 0.4, Pi/15];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 1:2 external resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{2,1}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{2,1}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{2,1}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{2,1}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 2:3 external resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{3,2}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{3,2}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{3,2}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{3,2}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 3:4 external resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{4,3}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{4,3}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{4,3}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{4,3}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 4:5 external resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{5,4}, 1.0, 0.1, Pi/20];
p2 = PlotPath[{5,4}, 1.0, 0.2, Pi/20];
p3 = PlotPath[{5,4}, 1.0, 0.3, Pi/20];
p4 = PlotPath[{5,4}, 1.0, 0.4, Pi/20];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]

Plot points along the paths of orbits at the 5:6 external resonance for eccentricities of 0.1, 0.2, 0.3 and 0.4.

p1 = PlotPath[{6,5}, 1.0, 0.1, Pi/15];
p2 = PlotPath[{6,5}, 1.0, 0.2, Pi/15];
p3 = PlotPath[{6,5}, 1.0, 0.3, Pi/15];
p4 = PlotPath[{6,5}, 1.0, 0.4, Pi/15];

Gather these using GraphicsArray.

Show[GraphicsArray[{p1,p2,p3,p4}]]