Chapter 8, Section 8.15, Figure 8.29a

First we input the masses (in kg) of the saturnian satellites Mimas, Enceladus, Tethys, Dione and Rhea using the data from Table A.9.

mass  =  {0.385,  0.73, 6.22, 10.52, 23.1}*10^(20)

Then we input their current semi-major axes again using data from Table A.9.

ao = {185520,  238020, 294660, 377400, 527040}

Defiine a function compare that calculates the ratio given in Eq.(8.251).  If this ratio is greater than 1 then the pair of satellites is approaching one another; if it is less than 1 then they are receding from one another..

compare[i_,  
        j_]  :=  (mass[[i]]/
            mass[[j]])*(ao[[j]]/ao[[i]])^(11/
              2)

Calculate and display the values of compare for all possible pairs of satellites.

Table[compare[i, j], {i, 1, 5}, {j, 1, 5}] // MatrixForm

Define a function a that returns the semi-major axis of satellite i at time t.  This is derived from Eq.(8.252) with K = 18.

a[i_, t_] := 
      N[Log[10, (ao[[i]]^(13/2) - (13/2)*18*mass[[i]]*
                  t)^(2/13)]]

Set the value of tmax, the maximum time in the past that we are considering.

tmax = 3.5*10^(12);

Derive a list of values of t and a for each satellite.

evol1  =  
        Table[{(tmax - t)/tmax, a[1, t]}, {t, 0, tmax, 
            0.01*10^(12)}];

evol2  =  
        Table[{(tmax - t)/tmax, a[2, t]}, {t, 0, tmax, 
            0.01*10^(12)}];

evol3  =  
        Table[{(tmax - t)/tmax, a[3, t]}, {t, 0, tmax, 
            0.01*10^(12)}];

evol4  =  
        Table[{(tmax - t)/tmax, a[4, t]}, {t, 0, tmax, 
            0.01*10^(12)}];

evol5  =  
        Table[{(tmax - t)/tmax, a[5, t]}, {t, 0, tmax, 
            0.01*10^(12)}];

Derive a plot of these values for each satellite.

p1 = 
      ListPlot[evol1, PlotJoined -> True,  
        DisplayFunction -> Identity]

p2 = 
      ListPlot[evol2,  PlotJoined -> True,  
        DisplayFunction -> Identity]

p3 = 
      ListPlot[evol3,  PlotJoined -> True,  
        DisplayFunction -> Identity]

p4 = 
      ListPlot[evol4,  PlotJoined -> True,  
        DisplayFunction -> Identity]

p5 = 
      ListPlot[evol5,  PlotJoined -> True,  
        DisplayFunction -> Identity]

Show all these plots together.

bigplot = 
      Show[p1, p2, p3, p4, p5, Frame -> True, 
        PlotRange -> {{0, 1}, {5.1, 5.8}}, 
        DisplayFunction -> $DisplayFunction,  AspectRatio -> 0.8]

Now we need to find the times when resonances have been encountered.  Define a function FindRes that will take the list of t and a values for two satellites and look for occurences of the ratio ratio.

\(FindRes[d1_, d2_, \ ratio_] := 
	
      Module[{}, 
		time\  = \ \(Transpose[
              d1]\)[\([1]\)]; 
		data1\  = \ 
          10^\((\(Transpose[d1]\)[\([2]\)])\); 
		data2\  = \ 
          10^\((\(Transpose[d2]\)[\([2]\)])\); 
		mm\  = \ 
          Abs[\((data1/data2)\)^\((3/2)\) - ratio]; 
		mindif\  = \ 
          Min[mm]; 
		tres\  = \ \(Flatten[
              Select[Transpose[{time, 
                    mm}], \ \((#[\([2]\)] == 
                      mindif)\) &]]\)[\([1]\)]; 
		pos\  = \ \(Flatten[
              Position[time, \ tres]]\)[\([1]\)]; 
		{d1[\([pos]\)], 
          d2[\([pos]\)]}]\)

Search for specific first- and second-order resonances.

l1 = FindRes[evol1, evol2,  N[1/2]]

l2 = FindRes[evol1, evol2,  N[2/3]]

l3  = FindRes[evol1, evol4, N[1/3]]

l4 = FindRes[evol2, evol3, N[5/6]]

l5 = FindRes[evol2, evol3, N[4/5]]

l6 = FindRes[evol2, evol3, N[3/4]]

l7 = FindRes[evol2, evol4, N[1/2]]

l8 = FindRes[evol3, evol4, N[2/3]]

l9 = FindRes[evol4, evol5, N[3/5]]

Create a graphic object consisting of lines for each of these resonances.

lines = 
      Graphics[Map[Line, {l1, l2, l3, l4, l5, l6, l7, l8, l9}]]

Display the final plot (see Fig.8.29a).

bigplot = 
      Show[p1, p2, p3, p4, p5, lines, Frame -> True, 
        PlotRange -> {{0, 1}, {5.1, 5.8}}, 
        DisplayFunction -> $DisplayFunction,  AspectRatio -> 0.8]