Chapter 10, Section 10.4, Figure 10.11

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

Define the resonance (2:1) and the eccentricity (0.1) for the first plot.

p=1/2; e=0.1;

Generate plots for 11 orbits in the rotating frame with the same eccentricity, linearly increasing semi-major axis and linearly decreasing (from zero) longitude of pericentre.

Do[
  ϖ = N[-i*Pi/12];
  a = 1.0 + i*0.15;
  rotlist={};
  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,4Pi,Pi/20}];
  plot[i] = 
  ListPlot[rotlist,AspectRatio->Automatic, 
         PlotJoined->True,Axes->None,
         DisplayFunction->Identity],
{i,0,10,1}]

Create graphics objects for the planet, its centre and the location of the exact resonance.

planet = Graphics[Circle[{0,0},0.4]];
resonance=
 Graphics[{Dashing[{0.02,0.02}],Circle[{0,0},0.8]}];
centre = Graphics[Point[{0,0}]];

Display the orbits, planet and resonance.  Note that for low eccentricity the orbits in the rotating frame for the 2:1 resonance (the m = 2 case) are centred ellipses.  (See Fig.10.11a.)

Show[planet, resonance, centre,
    Table[plot[i],{i,0,10,1}],
     DisplayFunction->$DisplayFunction,
     AspectRatio->Automatic]

Define the resonance (7:6) and the eccentricity (0.02) for the second plot.

p=6/7; e=0.02;

Generate plots for 21 orbits in the rotating frame with the same eccentricity, linearly increasing semi-major axis and linearly decreasing (from zero) longitude of pericentre.

Do[
  cp = N[-i*Pi/30];
  a = 1.0 + i*0.075;
  rotlist={};
  Do[
    posvec=position[t,a,e,cp];
    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,14Pi,Pi/10}];
  plot[i] = 
  ListPlot[rotlist,AspectRatio->Automatic, 
         PlotJoined->True,Axes->None,
         DisplayFunction->Identity],
{i,0,20,1}]

Display the orbits, planet and resonance.  Note that for low eccentricity the orbits in the rotating frame for the 7:6 resonance (the m = 7 case) produce seven-lobed shapes.  (See Fig.10.11b.)

Show[planet, resonance, centre,
 Table[plot[i],{i,0,20,1}],
     DisplayFunction->$DisplayFunction,
     AspectRatio->Automatic]