Chapter 9, Section 9.6, Figure 9.23

We start by defining the function f given in Eqn.(9.91).

f[ε_] := 1/Abs[(1+ε)^(-3/2)-1];

Now define the function ConjunctionMap that uses Eqs.(9.92)–(9.94) to take the system from one set of values of z, ε and λ to the next.

ConjunctionMap[mass_, g1_, {z_, ε_, λ_}] :=
Module[{znew, εnew, λnew},
  znew   = z + I*g1*Exp[I*λ]*mass;
  εnew = ε*Sqrt[1+4*(Abs[znew]^2 - Abs[z]^2)/(3*ε^2)];
  λnew = N[Mod[λ + 2*Pi*f[εnew], 2*Pi]];
  {znew, εnew, λnew}
]

Define a function, data, that will iterate the map n times returning the results in a list.

data[n_, mass_, {z_, ε_, λ_}] :=
Module[
{list = {{z, ε, λ}}, old = {z, ε, λ}, 
   ε1 = ε, 
   g1 = N[(8/9)*(2*BesselK[0,2/3]+
          BesselK[1,2/3])]*Sign[ε]/ε/ε},
Do[new = ConjunctionMap[mass, g1, old];
   old = new;
   AppendTo[list, new], 
  {n}];
list];

Define a function RegularChaotic that uses ConjunctionMap and the algorithm by Duncan et al. (1989) to determine if an orbit is regular (returning the value 0) or chaotic (returning the value 1).

RegularChaotic[mass_, ε_] :=
Module[
{old = {0, ε, 0}, ε1 = ε, e = 0, regcha = 0,
 g1 = N[(8/9)*(2*BesselK[0,2/3]+
        BesselK[1,2/3])]*Sign[ε]/ε/ε},
Do[new = ConjunctionMap[mass, g1, old];
   e = Abs[new[[1]]];
   If[e > ε1, regcha = 1];
   old = new, 
  {1000}];
regcha];

Define a function MaxEps that determines the maximum value of ε.

MaxEps[mass_] :=
Module[{},
upper = 2*1.24*mass^(2/7); 
lower = 0.5*1.24*mass^(2/7);
While[
   Abs[upper-lower]/lower >  0.01,
   middle = (upper + lower)/2;
   If[RegularChaotic[mass, middle] == 1,
      lower = middle,
      upper = middle]];
 middle];

Generate a list, mdata, of mass ratios that we will use.

mdata = Table[N[10^i], {i,-9,-3,0.1}];

Now calculate a list, edata, of the resulting values of ε.  N.B. This can take a very long time!

edata = Map[MaxEps[#] &, mdata];

Caculate logs of each list.

logmdata = Log[10, mdata];
logedata = Log[10, edata];

Create a list, logmedata, consisting of pairs of values of log mass and log ε

logmedata = Transpose[{logmdata, logedata}];

Plot this list.

pdata = 
ListPlot[logmedata, 
 PlotJoined->False,PlotRange->{{-9,-3},{-2.5,-0.5}},
 Frame->True]

Use the Mathematica Fit function to provide a linear fit to the data.

Fit[logmedata, {1, x}, x]

Define this fit as a function, fun0.

fun0[x_] := Fit[logmedata, {1, x}, x]

Define a function, fun1, to represent the relation given in Eq.(9.148).

fun1[x_] := Log[10,1.3] + N[2/7]*x

Plot the function fun0 ...

p0 = Plot[fun0[x],{x,-9,-3},
      PlotRange->{{-9,-3},{-2.5,-0.5}}, Frame->True]

... and the function fun1 ...

p1 = Plot[fun1[x],{x,-9,-3},
      PlotRange->{{-9,-3},{-2.5,-0.5}}, Frame->True]

... and put them on the same plot as the data points.  Note that the points and the fit may not be identical to those shown in Fig.9.23 because of the chaotic nature of the orbits.

Show[pdata,p0,p1]