Chapter 3, Section 3.6, Figure 3.7

First we define the function Location which returns the coordinates of the equlibrium 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}]

p1 = Location[0.2, 1]

p2 = Location[0.2, 2]

p3 = Location[0.2, 3]

p4 = Location[0.2, 4]

p5 = Location[0.2, 5]

Now put these positions together in a list and map the Point function on to each location.

lagpts=Map[Point,{p1,p2,p3,p4,p5}]

Now define a function czv that returns the value of the Jacobi constant of the zero velocity curve at the point (x,y).  (See Eq.(3.39).)

czv[x_,y_,μ_]:= 
2*((x^2 + y^2)/2 + μ/((-1 + μ + x)^2 + y^2)^(1/2) + 
   (1 - μ)/((μ + x)^2 + y^2)^(1/2));

We are now in a position to calculate the values of the Jacobi constant at the three critical points corresponding to the L1, L2 and L3.  We call these c1, c2 and c3.

c1 = czv[p1[[1]], p1[[2]], 0.2]

c2 = czv[p2[[1]], p2[[2]], 0.2]

c3 = czv[p3[[1]], p3[[2]], 0.2]

We then use ContourPlot to derive the zero-velocity curves for these three values of the Jacobi constant.

zvcon=
ContourPlot[czv[x,y,0.2],{x,-1.75,1.75},{y,-1.75,1.75},
            AspectRatio->Automatic,
            ContourShading->False,PlotPoints->200,
            Contours->{c1,c2,c3},
            Axes->False,DisplayFunction->Identity];

Now we can plot zvcon together with the locations of the five equilibrium points and points denoting the position of the two masses.

Show[zvcon,
     Graphics[{PointSize[0.02],lagpts,
              {Point[{-0.2,0}],Point[{0.8,0}]}}],
              Axes->False, Frame->False,
              DisplayFunction->$DisplayFunction]