Chapter 10, Section 10.3.1, Figure 10.7

The functions defined here are those already used in, for example, Fig.8.1–8.3.nb.

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]);

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_]:=
Module[{eca, x, y},
  eca = ea /. ksolve[ma,e];
  x = N[a(Cos[eca]-e)]; 
  y = N[a Sqrt[1-e^2] Sin[eca]];
  {x,y}]

We set the list of points, rotlist, to be empty.

rotlist={};

Define the resonance (3:2), the semi-major axis (unity) and the eccentricity (0.25).

p=3/2; a=1.0; e=0.25;

Create the list of points that make up the orbit.

Do[
(posvec=position[t,a,e];
f=N[ArcTan[posvec[[1]],posvec[[2]]]];
rad=N[Sqrt[posvec[[1]]^2+posvec[[2]]^2]];
th=p*t;
psi=th-f;
blip=N[{rad*Cos[psi],rad*Sin[psi]}];
AppendTo[rotlist,blip];),
{t,0,4Pi,Pi/20}]

Plot the points as the path in the rotating frame.

pa=ListPlot[rotlist,AspectRatio->Automatic, 
            PlotJoined->True,Axes->None]

Plot the points as individual points aloong the path in the rotating frame.

pb=
ListPlot[rotlist,AspectRatio->Automatic, Axes->None]

Show both plots together (see Fig.10.7).

Show[pa,pb,Prolog->PointSize[0.02]]