Chapter 8, Section 8.2, Figures 8.1–8.3

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

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

Now we define a function ksolve that solves Kepler's equation to find the eccentric 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 the function Orbit that plots an orbit of eccentricity, e, with semi-major axis specified by the rational made up from the list of two numbers given by pl.  The orbit of the outer planet is taken to be circular with unit semi-major axis.

Orbit[pl_List, e_] :=
Module[{time, p1, p2, p = Apply[Rational,pl],
        pt = If[pl[[1]] > pl[[2]], pl[[1]], pl[[2]]],
        sun=Graphics[{PointSize[0.1],Point[{0,0}]}],
        porb = Graphics[{Circle[{0,0},1]}],
        a = N[(pl[[1]]/pl[[2]])^(2/3)],
        circle = Graphics[{Circle[{0,0},1]}],
        box = 
        Graphics[Line[{{-1.1,-1.1},{1.1,-1.1},
                       {1.1,1.1},{-1.1,1.1},{-1.1,-1.1}}]]},
  ptime = Table[N[t], {t,0,pt*2*Pi,Pi/40}];
  posvec = Map[position[#,a,e]&,ptime];
  f = Map[N[ArcTan[#[[1]],#[[2]]]]&,posvec];
  rad = Map[N[Sqrt[#[[1]]^2+#[[2]]^2]]&,posvec];
  orbpts = Transpose[{rad*Cos[f],rad*Sin[f]}];
  p3 = 
  ListPlot[orbpts,
           PlotJoined->True,
           AspectRatio->Automatic, Axes->None, 
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  Show[box,sun,circle,p3,
       DisplayFunction->$DisplayFunction,
  AspectRatio->Automatic,PlotRange->All]]

Now plot the basic orbit for the 2:1 resonance used in Figs.8.1 and 8.2.

basic = Orbit[{1,2},0.5]

Duplicate this plot to form the figure.

Show[GraphicsArray[{basic,basic,basic,basic,basic}]]

Define a function PlotPathInertial to show the paths in the inertial frame around the unit orbit and the resonant orbit, adding points at equal time intervals, dt.

PlotPathInertial[pl_List, 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}]}],
        planet=Graphics[{PointSize[0.02],Point[{1.0,0.0}]}],
        porb = Graphics[{Circle[{0,0},1]}],
        a = N[(pl[[1]]/pl[[2]])^(2/3)],
        circle = Graphics[{Circle[{0,0},1]}],
        box = 
        Graphics[Line[{{-1.1,-1.1},{1.1,-1.1},
                       {1.1,1.1},{-1.1,1.1},{-1.1,-1.1}}]]},
  ptime = Table[N[t], {t,0,pt*2*Pi,Pi/40}];
  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];
  points2 = Transpose[{Cos[p*time],Sin[p*time]}];
  points1 = Transpose[{rad*Cos[f],rad*Sin[f]}];
  p1 = 
  ListPlot[points1,
           AspectRatio->Automatic, Axes->None,
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  p2 = 
  ListPlot[points2,
           AspectRatio->Automatic, Axes->None,
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  posvec = Map[position[#,a,e]&,ptime];
  f = Map[N[ArcTan[#[[1]],#[[2]]]]&,posvec];
  rad = Map[N[Sqrt[#[[1]]^2+#[[2]]^2]]&,posvec];
  orbpts = Transpose[{rad*Cos[f],rad*Sin[f]}];
  p3 = 
  ListPlot[orbpts,
           PlotJoined->True,
           AspectRatio->Automatic, Axes->None,
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  Show[box,sun,planet,circle,p1,p2,p3,
       DisplayFunction->$DisplayFunction,
  AspectRatio->Automatic,PlotRange->All]]

Plot the paths in the inertial frame for the 2:1 resonance with points at time intervals of π/4.  This gives the basic plot for Fig.8.3a.

p = PlotPathInertial[{1,2}, 0.2, Pi/4]

Define a function PlotPathRotate to show the paths in the rotating frame around the unit orbit and the resonant orbit, adding points at equal time intervals, dt.

PlotPathRotate[pl_List, e_, dt_] :=
Module[{p = Apply[Rational,Reverse[pl]],
        pt = If[pl[[1]] > pl[[2]], pl[[1]], pl[[2]]],
        sun=Graphics[{PointSize[0.05],Point[{0,0}]}],
        planet=Graphics[{PointSize[0.02],Point[{1.0,0.0}]}],
        porb = Graphics[{Circle[{0,0},1]}],
        a = N[(pl[[1]]/pl[[2]])^(2/3)],
        circle = Graphics[{Circle[{0,0},1]}],
        box = 
        Graphics[Line[{{-1.1,-1.1},{1.1,-1.1},
                       {1.1,1.1},{-1.1,1.1},{-1.1,-1.1}}]]},
  time = Table[N[t], {t,0,pt*2*Pi,dt}];
  ptime = Table[N[t], {t,0,pt*2*Pi,Pi/40}];
  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[f-time/p],rad*Sin[f-time/p]}];
  p1 = 
  ListPlot[points,
           AspectRatio->Automatic, Axes->None,
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  posvec = Map[position[#,a,e]&,ptime];
  f = Map[N[ArcTan[#[[1]],#[[2]]]]&,posvec];
  rad = Map[N[Sqrt[#[[1]]^2+#[[2]]^2]]&,posvec];
  orbpts = Transpose[{rad*Cos[f-ptime/p],rad*Sin[f-ptime/p]}];
  p3 = 
  ListPlot[orbpts,
           PlotJoined->True,
           AspectRatio->Automatic, Axes->None,
           DisplayFunction->Identity,
           PlotRange->{{-1.1,1.1},{-1.1,1.1}}];
  Show[box,sun,planet,p1,p3,
       DisplayFunction->$DisplayFunction,
  AspectRatio->Automatic,PlotRange->All]]

Plot the paths in the rotating frame for the 2:1 resonance with points at time intervals of π/4.  This gives the basic plot for Fig.8.3b.

p = PlotPathRotate[{1,2}, 0.2, Pi/4]