Chapter 3, Section 3.6, Figure 3.10

Again we make use of the function Location to allow us to calculate the position of Lagrangian point ptid for mass ratio μ.

Location[μ_, ptid_Integer] := 
      Module[{α, β, r1, r2, μ1 = 1 - μ, μ2 = μ, 
          lageq}, 
If[ptid == 1, 
          y = 0; 
α = (μ/(3*(1 - μ)))^(1/3); 
r20 = α - (1/
                    3)*α^2 - (1/
                    9)*α^3; 
lageq = μ1*(1 - 
                    r2 - 1/(1 - r2)^2) - μ2*(r2 - 
                    1/r2^2); 
r2 = 
            r2 /. FindRoot[
                lageq == 0, {r2, r20}]; 
r1 = 
            1 - r2; 
x = μ1 - 
              r2]; 
If[ptid == 2, 
          y = 0; 
α = (μ/(3*(1 - μ)))^(1/3); 
r20 = α + (1/
                    3)*α^2 - (1/
                    9)*α^3; 
lageq = μ1*(1 + 
                    r2 - 1/(1 + r2)^2) + μ2*(r2 - 
                    1/r2^2); 
r2 = 
            r2 /. FindRoot[
                lageq == 0, {r2, r20}]; 
r1 = 
            1 + r2; 
x = μ1 + 
              r2]; 
If[ptid == 3, 
          y = 0; 
β = -(7/
                      12)*(μ2/μ1) + (7/
                    12)*(μ2/μ1)^2; 
r10 = 
            1 + β; 
lageq = μ1*(r1 - 
                    1/r1^2) + μ2*(1 + r1 - 
                    1/(1 + r1)^2); 
r1 = 
            r1 /. FindRoot[
                lageq == 0, {r1, r10}]; 
r2 = 
            1 + r1;  x = -r1 - μ2]; 
If[
          ptid == 4, 
          x = 1/2 - μ2; 
y = 
            N[Sqrt[3]/2]; 
r1 = 1; 
r2 = 
            1]; 
If[ptid == 5, 
          x = 1/2 - μ2; 
y = 
            N[-Sqrt[3]/2]; 
r1 = 
            1; 
r2 = 1]; 
{x, 
          y}]

Create a list of logs of masses using the Table function.

logmass = Table[i, {i, -1, -10, -0.1}];

Now calculate the masses themselves for use in Location.

mass = 10^logmass;

Set up a loop to derive a ListPlot for each Lagrangian point where we plot the radial distance from the origin as a function of the mass.  Note we use the negative of logmass to get the mass decreasing from left to right.

Do[
p[i] = 
        ListPlot[
          Transpose[{-logmass, (Map[Location[#, i] &, 
                    mass])  /.  {xx_, yy_} ->  Sqrt[xx^2 + yy^2]}], 
          PlotJoined -> True, 
          DisplayFunction -> Identity], 
{i, 1, 
        5}]

Finally we display all the plots.  Note that the plots for L4 and L5 are identical.

Show[p[1], p[2], p[3], p[4], p[5], 
      PlotRange -> {{1, 10}, {0.9, 1.1}}, 
      DisplayFunction -> $DisplayFunction, AxesOrigin -> {1, 0.9}, 
      Frame -> True]