Chapter 3, Section 3.6, Figure 3.8

This is similar to the derivation of Fig.3.7 except that we are now illustrating the location of the points in the three-dimensional surface where the height is the negative of the Jacobi constant.

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}]

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));

For each equilibrium point we calculate the location of the point and the negative of the Jacobi constant of the zero velocity curve associated with that point, joining them together to create a point in three dimensions.

p1 = Location[0.2, 1];
pc1 = Join[p1, {-czv[p1[[1]], p1[[2]], 0.2]}]

p2 = Location[0.2, 2];
pc2 = Join[p2, {-czv[p2[[1]], p2[[2]], 0.2]}]

p3 = Location[0.2, 3];
pc3 = Join[p3, {-czv[p3[[1]], p3[[2]], 0.2]}]

p4 = Location[0.2, 4];
pc4 = Join[p4, {-czv[p4[[1]], p4[[2]], 0.2]}]

p5 = Location[0.2, 5];
pc5 = Join[p5, {-czv[p5[[1]], p5[[2]], 0.2]}]

Now we create a list of points from each of the five equilibrium points to use in the surface plot.

lagpts=Map[Point,{pc1,pc2,pc3,pc4,pc5}];

We derive the surface plot of the negative of the Jacobi constant for the zero velocity curves.

zvsurf=Plot3D[-czv[x,y,0.2],{x,-1.7,1.7},{y,-1.7,1.7},
       PlotPoints->50,Boxed->False,
       PlotRange->{-3.85,-2.8},
       ClipFill->None,Axes->None,Mesh->False,
       DisplayFunction->Identity];

Finally we add the points denoting the locations of the equilibrium points and display together with the surface plot.

Show[zvsurf,Graphics3D[{PointSize[0.02],lagpts}],
     DisplayFunction->$DisplayFunction,
     ViewPoint->{3.814, -1.510, 2.859}]