Chapter 3, Section 3.7, Figure 3.11

We start by defining 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)).  This was used in Fig.3.7.nb.

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

Check by finding the value of the Jacobi constant at the L4 point.

czv[1/2-0.1, Sqrt[3]/2, 0.1]

Produce a contour plot of the zero velocity curves for two values of the Jacobi constant above the L4 value.

cp = 
ContourPlot[czv[x,y,0.1],{x,-1.3,1.3},{y,0.2,1.5},
            ContourShading->False,Contours->{2.95,3.09},
            PlotPoints->100,AspectRatio->Automatic]

Display this together with a point to indicate the location of the L4 point.

Show[cp, Graphics[Point[{1/2-0.1,Sqrt[3]/2}]]]

We start by considering the larger curve.  Construct a table of x values that span the inner part of the curve.

data1i = Table[xx, {xx, -0.9, 0.9, 0.1}];

Now provide approximate y values for use as starting values in FindRoot later on.

start1i = 
{0.4, 0.513, 0.6, 0.66, 0.711, 0.748, 0.765, 0.78, 0.774,
 0.761, 0.736, 0.7, 0.65, 0.589, 0.517, 0.454, 0.42, 
 0.425, 0.467};

Use the same x values for the outer part of the curve and just add 0.4 (a typical width) to get estimates of the y values.

data1o = data1i;
start1o = start1i + 0.4;

Set up the left-hand side of the equation for each x that we need to solve.

test1 = Map[(Expand[czv[#,y,0.1]]-3.09)&, data1i];

Now solve the equations to get the set of points for the inner curve ...

inner1 = 
Table[{data1i[[i]], (y /. FindRoot[test1[[i]]==0,
           {y,start1i[[i]]}])},
       {i,1,Length[data1i]}];

... and the outer curve.

outer1 = 
Table[{data1o[[i]], (y /. FindRoot[test1[[i]]==0,
           {y,start1o[[i]]}])},
       {i,1,Length[data1o]}];

Add two extra points at the extremes for completeness.

extra11 = {-0.968, (y /. FindRoot[czv[-0.968,y,0.1]-3.09,
                        {y, 0.527}])}

extra12 = {0.953, (y /. FindRoot[czv[0.953,y,0.1]-3.09,
                        {y, 0.622}])}

Set up the inner and outer points to be plotted.

ipoints1 = Map[Point,inner1];

opoints1 = Map[Point,outer1];

Display these points plus the two extras on the main plot.

Show[cp,Graphics[{ipoints1,opoints1,Point[extra11],
Point[extra12]}]]

Define a function vector that will calculate a vector along the direction defined by (Overscript[x, ¨]) and (Overscript[y, ¨]) derived from the equations of motion.  Because \(TraditionalForm\`\(x\& . \)\)=\(TraditionalForm\`\(y\& . \)\)=0 on the zero velocity curve the initial velocity must be in the direction of ((Overscript[x, ¨]),(Overscript[y, ¨])).

vector[{x_,y_},μ_] :=
Module[{r1 = Sqrt[(x+μ)^2+y^2], r2 = Sqrt[(x-1+μ)^2+y^2]},
xvec = x - (1-μ)*(x+μ)/r1^3 - μ*(x-1+μ)/r2^3;
yvec = y - ((1-μ)/r1^3 + μ/r2^3)*y;
Line[{{x,y},{x,y}+{xvec, yvec}/5}]]

Calculate the vectors for all the points on the outer curve ...

vinner1 = Map[vector[#,0.1]&, inner1];
vouter1 = Map[vector[#,0.1]&, outer1];
vextra11 = vector[extra11,0.1];
vextra12 = vector[extra12,0.1];

... and plot them.

Show[cp,
  Graphics[{vextra11,vextra12,vinner1,vouter1}]]

Now we repeat the same procedure for the inner curve.  First we define a set of x values and give estimates of the corresponding y values for the inner and outer parts of the curve.

data2i = Table[xx, {xx, -0.05, 0.65, 0.1}];
start2i = 
{0.932, 0.89, 0.857, 0.82, 0.759, 0.713, 0.657, 0.643};
data2o = data2i;
start2o = {1.03, 1.05, 1.06, 1.04, 1.02, 0.97, 0.92, 0.829};

Set up the equations and solve them for the inner and outer parts of the curve.

test2 = Map[(Expand[czv[#,y,0.1]]-2.95)&, data2i];
inner2 = 
Table[{data2i[[i]], (y /. FindRoot[test2[[i]]==0,
           {y,start2i[[i]]}])},
       {i,1,Length[data2i]}];
outer2 = 
Table[{data2o[[i]], (y /. FindRoot[test2[[i]]==0,
           {y,start2o[[i]]}])},
       {i,1,Length[data2o]}];

Add two extra points for completeness.

extra21 = {-0.1, (y /. FindRoot[czv[-0.1,y,0.1]-2.95,
                        {y, 0.979}])}

extra22 = {0.7, (y /. FindRoot[czv[0.7,y,0.1]-2.95,
                        {y, 0.74}])}

Set up the vectors ...

vinner2 = Map[vector[#,0.1]&, inner2];
vouter2 = Map[vector[#,0.1]&, outer2];
vextra21 = vector[extra21,0.1];
vextra22 = vector[extra22,0.1];

... and plot them on the combined plot.

Show[cp,
  Graphics[{vextra11,vextra12,vinner1,vouter1,
            vextra21,vextra22,vinner2,vouter2}],
            AspectRatio->0.5]