Chapter 5, Section 5.4, Figure 5.13

Create a graphics object to denote the Sun.

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

We use the same functions, ksolve and position, as we used in Fig.5.9.nb, Fig.5.10a.nb, Fig.5.11.nb and Fig.5.12.nb to solve Kepler's equations and find the position vector as a function of mean anomaly.

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

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})

We will use a list, rotlist, to store the points we calculate.  First we must set it to be empty.

rotlist={};

We are interested in the p = 0 resonance with unit semi-major axis and an eccentricity, e = 0.5.

p=0; a=1.0; e=0.5;

Set up a loop to calculate and add points (positions at equal time intervals) to rotlist.  Note that we only have to go round the orbit once.

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

Plot the resulting points.

points=
ListPlot[rotlist,AspectRatio->Automatic, 
Prolog->PointSize[0.025],Axes->None]

Now go through the same procedure again but with more points in order to plot the path in the rotating frame.  First set the contents of rotlist to be empty.

rotlist = {};

Create the points in a loop ...

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

... and plot them.

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

Finally, display the points and the curve on the same plot.

Show[points,curve, sun]