Chapter 1, Section 1.7, Figure 1.8

We start by loading a standard Add-On Package that will be useful later on when we need to do continued fraction expansions of real numbers.

<<NumberTheory`ContinuedFractions`

Here we start the input of data on the various bodies to be included.  First we input the orbital periods of the planets in days (calculated using Kepler's third law) followed by their names.

planets = 365.25*
({0.38710353, 0.7233074, 1.00001161, 1.52371069,
  5.21021558, 9.53807012, 19.183302, 30.055144,
  39.53758}^(3/2));

planetsnames = {"Mercury", "Venus", "Earth", "Mars",
"Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"};

Now do the same for the jovian system ...

jupiter = 
 {0.49817905, 1.769137786, 3.551181041, 7.15455296, 
  16.6890184};

jupiternames = {"Amalthea", "Io", "Europa", "Ganymede",
"Callisto"};

... the saturnian system ...

saturn = 
 {0.942421813, 1.370217855, 1.887802160, 2.736914742,
  4.517500436, 15.94542068, 21.2766088, 79.3301825};

saturnnames = {"Mimas", "Enceladus", "Tethys", "Dione",
"Rhea", "Titan", "Hyperion", "Iapetus"};

... the uranian system ...

uranus = 
 {1.41347925, 2.52037935, 4.1441772, 8.7058717,
  13.4632389};

uranusnames = {"Miranda", "Ariel", "Umbriel",
"Titania", "Oberon"};

... and the neptunian system ...

neptune =
{5.8768541, 360.13619};

neptunenames = {"Triton", "Nereid"};

Set up a list containing the rationals that correspond to first-order commensurabilities of the form p/(p+1), with p≤100 and (1/3)added at the start.  This will form the set of rationals which will bound a ratio of orbital periods.

first = Join[{1/3},Table[i/(i+1), {i,1,100}]]

Define a function ratios which will create a list of ratios of periods >(1/3) + (1/12). The (1/12)is half of the "separation" between (1/3)and (1/2), the first two rationals in first.

\(ratios[data_List, \ name_List]\  := 

      Module[{len = Length[data], \ temp1, \ 
          temp2}, 
temp1 = 
          Table[{name[\([i]\)], name[\([j]\)], 
              data[\([i]\)]/data[\([j]\)]}, {i, 1, len - 1}, {j, i + 1, 
              len}]; 
temp2 = 
          FlattenAt[temp1, \ Table[{i}, {i, 1, len - 1}]]; 
        Select[temp2, \ \((#[\([3]\)] > 
                1/3 + 1/12)\) &]
]\)

Now create the relevant ratios for the planets and all the satellite systems.

prat = ratios[planets,  planetsnames]

jrat=ratios[jupiter, jupiternames]

srat=ratios[saturn, saturnnames]

urat=ratios[uranus, uranusnames]

nrat=ratios[neptune, neptunenames]

Define a function nearest that will find the two rationals from first that bound the each of the ratios of periods supllied.  nearest uses a continued fraction expansion of the ratio to find the closest rational from the list.  NB: The package NumberTheory`ContinuedFraction`  in Mathematica Version 4 is significantly different from that in Version 3.

nearest[vlist_List] :=
Module[{closest,pos,next},
closest = 
Map[
(Intersection[
Convergents[
ContinuedFractionForm[ContinuedFraction[#,4]]],
   first][[1]])&,
 vlist];
pos = Map[Position[first,#][[1,1]]&,closest];
next = 
Table[
  If[vlist[[i]] < closest[[i]], 
     first[[pos[[i]]-1]], 
     first[[pos[[i]]+1]]], {i,Length[vlist]}];
Map[Sort,Transpose[{closest,next}]]
]

Define a function ccalc to calculate the value of c (see Eqs.(1.16)Ð(1.19)), a measure of the normalised "distance" to the nearest first order resonance.

ccalc[data_List, names_List] :=
Module[{temp,rlist,nlist,upper,lower,a,b,c},
 temp = Transpose[ratios[data, names]];
 rlist = temp[[3]];
 nlist = nearest[rlist];
 upper = Transpose[nlist][[2]];
 lower = Transpose[nlist][[1]];
 a = (rlist-lower)/(upper-lower);
 b = Map[If[# <= 1/2, 0, 1]&, a];
 c = 2*Pi*(a-b);
 Transpose[{temp[[1]], temp[[2]], rlist, nlist, c}]
 ]

For example, here are the ratio of orbital periods, the bounding rationals and the values of c for valid pairs of planets.

ccalc[planets,  planetsnames]

Create a list of such results for all valid planet and satellite pairs.

clist = 
Union[ccalc[planets, planetsnames],ccalc[jupiter, jupiternames],ccalc[saturn, saturnnames],
      ccalc[uranus, uranusnames]]

Order this list by taking the absolute value of c.

orderedclist = 
      Sort[clist,  (Abs[#2[[5]]]  >  
              Abs[#1[[5]]]) &]

Isolate just the values of c.

c  =  Transpose[orderedclist][[5]]

Note the total number of valid ratios.

n  =  Length[c]

Create a base line going from 0 to π ...

baseline=Graphics[{Line[{{0,-0.125},{Pi,-0.125}}]}];

... and then create a line for each of the pairs.

lines=Graphics[Map[Line[{{#,0},{#, 0.25}}]&,Abs[c]]];

Now display the main plot ...

Show[baseline,lines,
AspectRatio->9/80,PlotRange->{{0,Pi},{-0.13,0.25}}]

... and the close-up of the region in the range 0 to 0.15.

Show[baseline,lines,
AspectRatio->9/80, PlotRange->{{0,0.15},{-0.13,0.25}}]