This notebook is about the planar, circular restricted three-body problem.  There are two primary masses P1 of size mu and P2 of size 1-mu  where 0<mu<1.  In rotating coordinates they are located at positions (1-mu,0) and (-mu,0) respectively.  A third mass moves under the gravitational influence of the two primaries.

Solving the differential equation and plotting some orbits

The differential equation are given in terms of an amended potential V which encompasses both the gravitational interactions and the centrifugal force due to the rotating coordinates.  Here is the potential and its partitial derivatives.

V = (x^2+y^2)/2+ (1-mu)/Sqrt[(x+mu)^2+y^2]+mu/Sqrt[(x+mu-1)^2+y^2]
Vx=D[V,x]
Vy = D[V,y]

Here are the differential equations .  The velocity variables are u, v.  The equations reflect gravity, the centrifugal force and the Coriolus force.

xdot = u[t];
ydot = v[t];
udot = 2*v[t]+Vx/.{x->x[t],y->y[t]};
vdot = -2*u[t]+Vy/.{x->x[t],y->y[t]};
r3bpvf = {xdot, ydot, udot, vdot}

Here is how to set up the equations for solving with Mathematica.

r3bpode ={x'[t]==xdot,y'[t]==ydot, u'[t]==udot,v'[t]==vdot}

The following function  solves the ODE with given initial conditions at time 0 and mass mu over the time interval t0, t1

MakeR3BPOrbit[{x0_,y0_,u0_,v0_},mu0_,t0_,t1_]:=
NDSolveValue[(r3bpode/.mu->mu0)~Join~{x[0]==x0,y[0]==y0,u[0]==u0,v[0]==v0},{x[t],y[t],u[t],v[t]},{t,t0,t1}]

Here is  a function to plot the (x,y) coordinates of the solutions in the rotating system.  The primaries are shown too.

MakeR3BPPlot[orb_,mu0_,t0_,t1_,options___]:=Module[{pp,p1,p2},
pp = ParametricPlot[orb[[{1,2}]],{t,t0,t1},options,PlotRange->{{-1,1},{-1,1}}];
p1= Graphics[{Blue,PointSize[0.03],Point[{-mu0,0}]}];
p2= Graphics[{Green,PointSize[0.03],Point[{1-mu0,0}]}];
Show[pp,p1,p2]]

Here is an example.

mu0=1/3.;
x0=0.5;y0=0;
u0=0;v0=1;
t10=10;
orb1= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t10]
MakeR3BPPlot[orb1,mu0,0,10]

This orbit is moving around the primary with mass mu.  To see the actual motion we can make a plot in the nonrotating, inertial coordinate system.  We need to plot three curves and three points.

MakeR3BPPlotInertial[orb_,mu0_,t0_,t1_,options___]:=Module[{pp,p0,p1,p2,c0,c1,c2},
c0 = {Cos[t]*orb[[1]]-Sin[t]*orb[[2]],Sin[t]*orb[[1]]+Cos[t]*orb[[2]]};
c1={Cos[t]*(-mu0),Sin[t]*(-mu0)};
c2={Cos[t]*(1-mu0),Sin[t]*(1-mu0)};
pp = ParametricPlot[{c0,c1,c2},{t,t0,t1},PlotRange->{{-1,1},{-1,1}},PlotStyle->{Black,Blue,Green,}options];
p0= Graphics[{Black,PointSize[0.03],Point[c0/.t->t1]}];p1= Graphics[{Blue,PointSize[0.03],Point[c1/.t->t1]}];
p2= Graphics[{Green,PointSize[0.03],Point[c2/.t->t1]}];
Show[pp,p0,p1,p2]]

MakeR3BPPlotInertial[orb1,mu0,0,t10]

To see an animation, use the Manipulate command.  You can compare the motion in the rotating and intertial frames.  You can see that for this orbit, the “moon” is orbiting the “planet” in a  retrograde fashion.  This example also introduces a way to make a “worm plot” where only the path for a short time is shown.

Manipulate[GraphicsGrid[{{MakeR3BPPlot[orb1,mu0,0,t1],MakeR3BPPlotInertial[orb1,mu0,Max[0.0001,t1-1],t1]}}],{t1,0.001,t10}]

Here is an interactive plotter where you can select the initial position and velocity and then see the orbit.

MakeR3BPOrbitAndPlot[{x0_,y0_,u0_,v0_},mu0_,t0_,t1_,options___]:= Module[{orb},
orb = MakeR3BPOrbit[{x0,y0,u0,v0},mu0,t0,t1];
Show[MakeR3BPPlot[orb,mu0,t0,t1,options],Graphics[Arrow[{{x0,y0},{x0+0.3*u0,y0+0.3*v0}}]]]]

MakeR3BPOrbitAndPlotInertial[{x0_,y0_,u0_,v0_},mu0_,t0_,t1_,options___]:= Module[{orb},
orb = MakeR3BPOrbit[{x0,y0,u0,v0},mu0,t0,t1];
Show[MakeR3BPPlotInertial[orb,mu0,t0,t1,options],Graphics[Arrow[{{x0,y0},{x0+0.3*u0,y0+0.3*v0}}]]]]

Manipulate[MakeR3BPOrbitAndPlot[{q[[1]],q[[2]],v[[1]],v[[2]]},mu0,0,tmax0,PlotRange->{{-size,size},{-size,size}}],{{tmax0,0.001},0.001,20,ImageSize->Large},{{mu0,0.5},0,1,ControlPlacement->Right,ImageSize->Small},
{{size,3},0.1,10,ControlPlacement->Right,ImageSize->Small},{{q,{1,0}},Locator},{{q,{1-mu0,0}},InputField,ControlPlacement->Right},{{v,{0,1}},InputField,ControlPlacement->Right},{{v,{0,1}},{-1,-1},{1,1},Appearance->"Labeled",ControlPlacement->Right,ImageSize->{150,150}}]

Manipulate[MakeR3BPOrbitAndPlotInertial[{q[[1]],q[[2]],v[[1]],v[[2]]},mu0,Max[0.0001,tmax0-1],tmax0,PlotRange->{{-size,size},{-size,size}}],{{tmax0,0.001},0.001,20,ImageSize->Large},{{mu0,0.5},0,1,ControlPlacement->Right,ImageSize->Small},
{{size,3},0.1,10,ControlPlacement->Right,ImageSize->Small},{{q,{1,0}},Locator},{{q,{1-mu0,0}},InputField,ControlPlacement->Right},{{v,{0,1}},InputField,ControlPlacement->Right},{{v,{0,1}},{-1,-1},{1,1},Appearance->"Labeled",ControlPlacement->Right,ImageSize->{150,150}}]

Hill’s regions and Lagrange points

For a given energy, the position (x,y) is confined to the HIll’s region where V +h >= 0.  It’s easy to plot these.

MakeHillPlot[mu0_,h0_,options___] :=
RegionPlot[(V/.mu->mu0)+h0≥ 0,{x,-2,2},{y,-2,2},PlotStyle->{LightBlue},PlotPoints->30,options]

MakeHillPlot[1/3,-1.9]

Manipulate[MakeHillPlot[mu0,h0],{{mu0,1/2},0,1},{{h0,-2},-4,0}]

To show an orbit insides its Hill’s region, you need to match the energy of the orbit to the Hill’s region.

Energy[{x0_,y0_,u0_,v0_},mu0_]:=(u0^2+v0^2)/2-V/.{x->x0,y->y0,mu->mu0}

Here is the orbit used as an example in above.

mu0=1/3.;
x0=0.5;y0=0;
u0=0;v0=1;
t10=10;
h0 = Energy[{x0,y0,u0,v0},mu0];
orb1= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t10];
hillplot1=MakeHillPlot[mu0,h0];
orbplot1= MakeR3BPPlot[orb1,mu0,0,t10];
Show[hillplot1,orbplot1,PlotRange->{{-1,1},{-1,1}}]

If you want to find an orbit in a certain Hill’s region, you have to choose an initial condition with the right energy.  Here is a  function which sets the energy by scaling the given velocity vector.  As long as the initial position is in the Hill’s region, it returns an initial condition with the right energy.

SetEnergy[{x0_,y0_,u0_,v0_},mu0_,h0_]:=Module[{vnormsq,u1,v1},
vnormsq = 2*(V+h0)/.{x->x0,y->y0,mu->mu0}; (*the desired velocity norm squared *)
If [vnormsq<0, Print["Position not in Hill's region"];Return[{x0,y0,u0,v0}]];
{u1,v1}={u0,v0}*Sqrt[Abs[vnormsq]]/Sqrt[u0^2+v0^2];
{x0,y0,u1,v1}]

Here’s an orbit in a higher energy level.  This  could be a spacecraft moving from near  the Earth to near the Moon.

mu0=1/3.;
x0=0.2;y0=0;
u0=0;v0=1;
t12=30;
h0 =-1.9;
{x0,y0,u0,v0}= SetEnergy[{x0,y0,u0,v0},mu0,h0]
orb2= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t12];
hillplot2=MakeHillPlot[mu0,h0];
orbplot2= MakeR3BPPlot[orb2,mu0,0,t12];
Show[hillplot2,orbplot2,PlotRange->{{-1,1},{-1,1}}]
Manipulate[GraphicsGrid[{{MakeR3BPPlot[orb2,mu0,0,t1],MakeR3BPPlotInertial[orb2,mu0,Max[0.0001,t1-1],t1]}}],{t1,0.001,t12}]

The zero velocity curves for a given mu and all energies are just the level curves of the potential V.

GraphicsGrid[{{Plot3D[(V/.mu->mu0),{x,-2,2},{y,-2,2}],ContourPlot[(V/.mu->mu0),{x,-2,2},{y,-2,2},ContourShading->False,Contours->30,PlotPoints->30]}}]

The critical points of V are the five Lagrange points.  L4 and L5 are the equilateral points.

L4={1/2-mu,Sqrt[3]/2};L5={1/2-mu,-Sqrt[3]/2};
L45energy = Energy[L4~Join~{0,0},mu]//Simplify

The others are collinear and their approximate positions can be found numerically for a given value of mu.  Here is the graph of the potential over the x axis for mu = 1/3.  You can see the three critical points.

V0=(V/.{y->0,mu->1/3})
Plot[V0,{x,-2,2}]

The function FindMinimum will give the position of the critical point and the value of the potential there.

mu0=1/3;
min1=FindMinimum[(V /. {y -> 0, mu -> mu0}), {x,.4}]
L1={x/.min1[[2]],0}
L1energy = -min1[[1]]

min2=FindMinimum[(V /. {y -> 0, mu -> mu0}), {x, -1}]
L2={x/.min2[[2]],0}
L2energy = -min2[[1]]

min3=FindMinimum[(V /. {y -> 0, mu -> mu0}), {x, 1}]
L3={x/.min3[[2]],0}
L3energy = -min3[[1]]

Using this information, you can plot the most interesting zero velocity curves and show the Lagrange points.

mu0=1/3;cp = ContourPlot[(V/.mu->mu0),{x,-1.6,1.6},{y,-1.6,1.6},ContourShading->False,Contours->{-L1energy,-L2energy,-L3energy,-1.1*L45energy/.mu->mu0,-2*L1energy/.mu->mu0},PlotPoints->30];
lagrangepoints = Graphics[{PointSize[0.03],Black, Point[{L1,L2,L3,L4/.mu->mu0,L5/.mu->mu0}]}];
Show[cp,lagrangepoints]

In rotating coordinates, these are all equilibrium points.  In nonrotating coordinates, they represent periodic orbits.  For mu = 1/3 they are all unstable.  To simulate this, you can just write down the Lagrange orbits explicitly.  For example, here are the orbit  at L1 and L4.  All three masses rotate at the same speed, maintaining a constant shape.

L1orbit = L1~Join~{0,0}
L4orbit = (L4/.mu->mu0)~Join~{0,0}
Manipulate[GraphicsGrid[{{MakeR3BPPlotInertial[L1orbit,mu0,Max[0.0001,t1-1],t1],MakeR3BPPlotInertial[L4orbit,mu0,Max[0.0001,t1-1],t1]}}],{t1,0.001,8*Pi}]

Stability of Lagrange points and Lyapunov orbits

To analyze the stability of the Lagrange points, we need to find the eigenvalues and eigenvectors of the linearized ODE.  Here is the ODE and its derivative matrix.

f=r3bpvf/.{x[t]->x,y[t]->y,u[t]->u,v[t]->v}
df = {D[f,x],D[f,y],D[f,u],D[f,v]}//Transpose

First consider the equilateral points.  By symmetry, L4 and L5 have the same eigenvalues.   As examples, consider the equal mass case and the case of the Earth-Moon system, mu = 1/82.  In the first case there are four nonreal eigenvalues  with nonzero real parts and the equilibrium is unstable.  In the second case, all eigenvalues are  imaginary and L4 is linearly stable.

dfL4= df/.{x->1/2-mu,y->Sqrt[3]/2,u->0,v->0}//Simplify;
dfL4//MatrixForm
CharacteristicPolynomial[dfL4,z]
evL4= Eigenvalues[dfL4]
evL4/.{mu->1/2}
evL4/.{mu->1/82}

Consider the linearly stable case for  mu = 1/82.   Although the eigenvalues do not imply true nonlinear stability, they do show us that the nearby orbits are approximately a linear combination of two oscillations.  Here is a movie of a typical nearby orbit starting a bit to the left of L4.  Perhaps there is one of the Trojan asteroids which behaves like this.

mu0=1/82.;
x0=1/2-mu0-0.01;y0=Sqrt[3]/2;
u0=0;v0=0;
t10=100;
h0 = Energy[{x0,y0,u0,v0},mu0];
orblyap= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t10];
Manipulate[Show[MakeR3BPPlot[orblyap,mu0,Max[0.0001,t1-10],t1],Graphics[{Black,PointSize[0.02],Point[{1/2-mu0,Sqrt[3]/2}]}],PlotRange->{{-.1,1.1},{-.1,1.1}}],{t1,0.00001,t10}]

Now for some periodic motions.  Each imaginary pair determines an invariant plane filled with periodic solutions of the linearized ODE.   According to the Lyapunov center theorem, there should really be two families of periodic orbits close to these planes.  This will be explored numerically for  mu = 1/82.  First find the invariant planes by finding the real and imaginary parts of the eigenvectors.  The Chop command eliminates spurious roundoff errors.

evL40 = evL4/.mu->1/82.
evectL4=Eigenvectors[dfL4/.mu->1/82.]//Simplify
z1=evectL4[[1]]//Chop;
x1=Re[z1]
y1=Im[z1]
z2=evectL4[[3]]//Chop;
x2=Re[z2]
y2=Im[z2]

By moving a small distance away from L4 inside one of the eigenplanes, we can find an approximation to the initial condition for the Lyapunov orbit.  While the resulting orbit is only approximately periodic, there is an exactly periodic orbit nearby.  How to find them is rather a long story which will be told later.

mu0=1/82.;
eps0 = 0.1;
{dx,dy,du,dv}=eps0*y1;
x0=1/2-mu0+dx;y0=Sqrt[3]/2+dy;
u0=du;v0=dv;
t10=25;
h0 = Energy[{x0,y0,u0,v0},mu0];
orblyap4= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t10];
hillplotlyapl4=MakeHillPlot[mu0,h0];
lyapplot4= MakeR3BPPlot[orblyap4,mu0,0,t10];
Show[hillplotlyapl4,lyapplot4,Graphics[{Black,PointSize[0.02],Point[{1/2-mu0,Sqrt[3]/2}]}],PlotRange->{{-.1,1.1},{-.1,1.1}}]

Instead, let’s have a look at the collinear Lagrange point L1 for mu = 1/3.  We have to find L1 numerically, compute the eigenvalues, etc.

Plot[V /. {y -> 0, mu -> mu0},{x,-1,2}]

mu0=1/3;
min1=FindMinimum[(V /. {y -> 0, mu -> mu0}), {x,.3}]
L1={x/.min1[[2]],0}
L1energy = -min1[[1]]

dfL1= df/.{x->L1[[1]],y->L1[[2]],u->0,v->0,mu->mu0}//Simplify;
dfL1//MatrixForm
CharacteristicPolynomial[dfL1,z]
evL1= Eigenvalues[dfL1]//Chop
evectL1=Eigenvectors[dfL1/.mu->1/3.]//Simplify
z3=evectL1[[3]]//Chop;
x3=Re[z3] (* eigenspace of the imaginary pair *)
y3=Im[z3]

Although L1 is unstable, there is still a Lyapunov family of periodic orbit nearby corresponding to the imaginary pair of eigenvalues. Here is an attempt to choose an initial condition in the eigenspace of the imaginary pair of eigenvalues.   The resulting orbit  does not look even approximately periodic.   In fact  it quickly moves away from the eigenspace due to the real, highly unstable eigenvalue  3.73009.

mu0=1/3.;
eps0 = 0.2;
{dx,dy,du,dv}=eps0*x3;
{x0,y0,u0,v0}= (L1~Join~{0,0})+{dx,dy,du,dv}
t10=5;
h0 = Energy[{x0,y0,u0,v0},mu0]
orblyap1= MakeR3BPOrbit[{x0,y0,u0,v0},mu0,0,t10];
hillplotlyapl1=MakeHillPlot[mu0,h0,PlotRange->{{0,0.5},{-.25,.25}},PlotPoints->30];
lyapplot1= MakeR3BPPlot[orblyap1,mu0,0,t10];
Show[hillplotlyapl1,lyapplot1,Graphics[{Black,PointSize[0.02],Point[L1]}],PlotRange->{{0,0.5},{-.25,.25}}]

Finding the Lyapunov orbits

To find a more accurate initial conditions for the Lyapunov orbits, we will try to solve for the initial condition as a fixed point of a Poincare section.  This section is rather complicated and you may want to just run it and look at the results for now.

Let’s find the Lyapunov families near L4.  Since the approximate initial  point eps0*y1 has dx = 0, it is convenient to use x = 1/2 - mu0 as a section.   Since du <0, orbits leave the section in the direction of decreasing x.  We will look for the next time the orbit reaches x = 1/2-mu0 with x decreasing.  This is done using  Mathematica’s EventLocator capability with the option “Direction”->-1 below.

MakeOrbitToL4Section[{x0_,y0_,u0_,v0_},mu0_,tmax_]:=NDSolveValue[(r3bpode/.mu->mu0)~Join~{x[0]==x0,y[0]==y0,u[0]==u0,v[0]==v0},{x[t],y[t],u[t],v[t]},{t,0,tmax},
Method->{ "EventLocator","Event"->x[t]-1/2+mu0,"Direction"->-1 }]

Next we use this to make a 3D Poincare map function and a 2D version with fixed energy.  The latter requires setting the energy by changing only v.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

PoincareMapL4Section[{y0_,u0_,v0_},mu0_,tmax_]:=Module[{orb, tfinal},
orb = MakeOrbitToL4Section[{1/2-mu0,y0,u0,v0},mu0,tmax];
tfinal =InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
orb[[{2,3,4}]]/.t->tfinal (*read off y,u,v at the final time *)
]

SetEnergy2[{x0_,y0_,u0_,v0_},mu0_,h0_]:={x0,y0,u0,Sign[v0]*Sqrt[2((V/.{x->x0,y->y0,mu->mu0})+h0)-u0^2]}

PoincareMapL4Section[{y0_,u0_},mu0_,h0_,tmax_]:=Module[{orb, tfinal,x1,y1,u1,v1},
{x1,y1,u1,v1} =  SetEnergy2[{1/2-mu0,y0,u0,1},mu0,h0];
orb = MakeOrbitToL4Section[{x1,y1,u1,v1},mu0,tmax];
tfinal =InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
orb[[{2,3}]]/.t->tfinal (*read off y,u at the final time *)
]

To use Newton’s method for finding fixed points, we need the Jacobian matrix of the Poincare map.  This can be approximated by just  evaluating the Poincare map at a couple  of very close points.  Here is the 2D, fixed energy version.

PoincareMapL4SectionWithJacobian[{y0_,u0_},mu0_,h0_,tmax_,delta_:0.000001]:=Module[{orb, tfinal,x1,y1,u1,v1,yf,uf,yff,uff,a,b,c,d},
{x1,y1,u1,v1} =  SetEnergy2[{1/2-mu0,y0,u0,1},mu0,h0];
orb = MakeOrbitToL4Section[{x1,y1,u1,v1},mu0,tmax];
tfinal =InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
{yf,uf} =orb[[{2,3}]]/.t->tfinal ; (* final point *)
{x1,y1,u1,v1} =  SetEnergy2[{1/2-mu0,y0+delta,u0,1},mu0,h0];
orb = MakeOrbitToL4Section[{x1,y1,u1,v1},mu0,tmax];
tfinal =InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
{yff,uff} =orb[[{2,3}]]/.t->tfinal ; (* final point *)
{a,c} = {yff-yf,uff-uf}/delta; (* approx to  y partial derivatives *)
{x1,y1,u1,v1} =  SetEnergy2[{1/2-mu0,y0,u0+delta,1},mu0,h0];
orb = MakeOrbitToL4Section[{x1,y1,u1,v1},mu0,tmax];
tfinal =InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
{yff,uff} =orb[[{2,3}]]/.t->tfinal ; (* final point *)
{b,d} = {yff-yf,uff-uf}/delta; (* approx to  u partial derivatives *)
{{yf,uf},{{a,b},{c,d}}}  (* image point and approximate derivative matrix *)
]

Now use Newton’s method to find the fixed point.  The second version returns the initial condition {0, y,u,v} of the fixed point {y,u}.

FindPoincareMapL4FixedPoint[{y0_,u0_},mu0_,h0_,tol_:10^(-10),maxiters_:100]:= Module[{y1,u1,y2,u2,a,b,c,d,dely,delu,n=0},
{y1,u1}= {y0,u0};
While[n<maxiters,
{{y2,u2},{{a,b},{c,d}}} = PoincareMapL4SectionWithJacobian[{y1,u1},mu0,h0,100];
dely = y2-y1;delu=u2-u1;
If[Sqrt[dely^2+delu^2]≤tol,Return[{y2,u2}]];
{y1,u1} = {y1,u1} - {(d-1)*dely-b*delu,-c*dely+(a-1)*delu}/((a-1)(d-1)-b*c);
n++;];
{y0,u0}]

FindPoincareMapL4FixedPoint2[{y0_,u0_},mu0_,h0_,tol_:10^(-10),maxiters_:100]:= Module[{y1,u1,y2,u2,a,b,c,d,dely,delu,n=0},
{y1,u1}= {y0,u0};
While[n<maxiters,
{{y2,u2},{{a,b},{c,d}}} = PoincareMapL4SectionWithJacobian[{y1,u1},mu0,h0,100];
dely = y2-y1;delu=u2-u1;
If[Sqrt[dely^2+delu^2]≤tol,Return[SetEnergy2[{1/2-mu0,y2,u2,1},mu0,h0]]];
{y1,u1} = {y1,u1} - {(d-1)*dely-b*delu,-c*dely+(a-1)*delu}/((a-1)(d-1)-b*c);
n++;];
{y0,u0}]

Here is a plot of the orbit corresponding to our fixed point.  The orbit moves clockwise around L4.

mu0=1/82.;
eps0 = 0.1;
{dx,dy,du,dv}=-eps0*y2;
x0=1/2-mu0+dx;y0=Sqrt[3]/2+dy;
u0=du;v0=dv;
{x0,y0,u0,v0}
h0 = Energy[{x0,y0,u0,v0},mu0]
PoincareMapL4Section[{y0,u0},mu0,h0,10]
{x00,y00,u00,v00}=FindPoincareMapL4FixedPoint2[{y0,u0},mu0,h0]
Show[Table[MakeR3BPPlot[MakeR3BPOrbit[i{x00,y00,u00,v00},mu0,0,30],mu0,0,30,PlotRange->{{-.3,1.1},{-.3,1.1}},PlotStyle->Blue],{i,1,10}]~Join~{Graphics[{Black,PointSize[0.02],Point[{1/2-mu0,Sqrt[3]/2}]}]}]

To look at the whole family of Lyapunov orbits, just repeat this process of finding the fixed point with several initial guesses.  The energies are a bit above that of L4 (less negative).

mu0=1/82.;
guesstable= Table[{1/2-mu0,Sqrt[3]/2,0,0}+eps*y1,{eps,0.01,0.501,0.05}];
energytable=Map[Energy[#,mu0]&,guesstable]
ictable=Table[FindPoincareMapL4FixedPoint2[guesstable[[i,{2,3}]],mu0, energytable[[i]]],{i,1,Length[guesstable]}];
l4lyapplot1=Show[Table[MakeR3BPPlot[MakeR3BPOrbit[ictable[[i]],mu0,0,30],mu0,0,30,PlotRange->{{-.3,1.1},{-.3,1.1}},PlotStyle->Blue],{i,1,10}]~Join~{Graphics[{Black,PointSize[0.02],Point[{1/2-mu0,Sqrt[3]/2}]}]}]

To find the other Lyapunov family, use the eigenspace of the other imaginary pair.  The energies here are below that of L4.

mu0=1/82.;
eps0 = 0.1;
guesstable= Table[{1/2-mu0,Sqrt[3]/2,0,0}-eps*y2,{eps,0.01,0.501,0.05}];
energytable=Map[Energy[#,mu0]&,guesstable]
ictable=Table[FindPoincareMapL4FixedPoint2[guesstable[[i,{2,3}]],mu0, energytable[[i]]],{i,1,Length[guesstable]}];
l4lyapplot2=Show[Table[MakeR3BPPlot[MakeR3BPOrbit[ictable[[i]],mu0,0,30],mu0,0,30,PlotRange->{{-.3,1.1},{-.3,1.1}},PlotStyle->Green],{i,1,10}]~Join~{Graphics[{Black,PointSize[0.02],Point[{1/2-mu0,Sqrt[3]/2}]}]}]

Here are both families.

Show[GraphicsGrid[{{l4lyapplot1,l4lyapplot2}}]]

Levi Civita regularization

To handle near collisions  with the  primary P2 of mass mu at position (1-mu,0) a translation of the complex squaring map can be used.  The velocity vector  also gets transformed.  New variables are xi, eta, a, b

phi[{xi_,eta_},mu_] :={ 1-mu +xi^2-eta^2,2*xi*eta}
phi[{xi_,eta_,a_,b_},mu_] :=Module[{x,y},
{x,y} = { 1-mu +xi^2-eta^2,2*xi*eta};
{x,y,2*(xi*a-eta*b)/(xi^2+eta^2),2*(xi*b+eta*a)/(xi^2+eta^2)}]

Since the tranformation is two-to-one, there are two inverse  branches.  These will  be based on  the branches of the complex square root.

ComplexSqrt1[{xi_,eta_}]:= Module[{r=Sqrt[xi^2+eta^2]},If[eta==0&&xi<0,{0,Sqrt[(r-xi)/2]},{Sqrt[(r+xi)/2],Sign[eta]*Sqrt[(r-xi)/2]}]]
ComplexSqrt2[{xi_,eta_}]:=-ComplexSqrt1[{xi,eta}]

phiinverse1[{x_,y_},mu_]:=ComplexSqrt1[{x-1+mu,y}]
phiinverse1[{x_,y_},mu_]:=ComplexSqrt2[{x-1+mu,y}]
phiinverse1[{x_,y_,u_,v_},mu_]:=Module[{xi,eta,a,b},
{xi,eta} = ComplexSqrt1[{x-1+mu,y}];
{a,b} = {(xi*u+eta*v)/2,(xi*v-eta*u)/2};
{xi,eta,a,b}]
phiinverse2[{x_,y_,u_,v_},mu_]:=Module[{xi,eta,a,b},
{xi,eta} = ComplexSqrt2[{x-1+mu,y}];
{a,b} = {(xi*u+eta*v)/2,(xi*v-eta*u)/2};
{xi,eta,a,b}]

Here  is the regularized potential W(xi,eta) and its partial derivatives.  Also the energy  in regularized variables.

WLeviCivita = (xi^2+eta^2)((xi^2-eta^2+1-mu)^2+4*xi^2*eta^2)/8 +(1-mu)*(xi^2+eta^2)/(4*Sqrt[(xi^2-eta^2+1)^2+4*xi^2*eta^2])+(mu+h*(xi^2+eta^2))/4
WLCxi=D[WLeviCivita,xi]
WLCeta = D[WLeviCivita,eta]

LeviCivitaEnergy[{xi0_,eta0_,a0_,b0_},mu0_]:=4*((a0^2+b0^2)/2-(WLeviCivita/.{xi->xi0,eta->eta0,mu->mu0,h->0}))/(xi0^2+eta0^2)

mu0=mu;
xi0=xi;eta0=0;
W0=WLeviCivita/.{xi->xi0,eta->eta0,mu->mu0,h->0}//Simplify[#,Assumptions->0<mu<1]&
{x0,y0}=phi[{xi0,eta0},mu0]
V0=V/.{x->x0,y->y0,mu->mu0}//Simplify[#,Assumptions->0<mu<1]&;
(xi0^2+eta0^2)*V0/4//Expand
W0//Expand
%-%%

mu0=1/3;
xi0=0.2;eta0=0;
a0=0;b0=1.;
{x0,y0,u0,v0}=phi[{xi0,eta0,a0,b0},mu0]
h0=LeviCivitaEnergy[{xi0,eta0,a0,b0},mu0]
Energy[{x0,y0,u0,v0},mu0]
phiinverse1[{x0,y0,u0,v0},mu0]

Here is the regularized vectorfield, ODE in Mathematica notation and a function  for making an orbit.  This time we need to calculate  the energy and substitute it into  the regularized ODE.  In addition to finding the regularized positions and velocities, the unregularized time will be found as the 5th component of the solution.

xidot = a[t];
etadot = b[t];
adot = 2*(xi^2+eta^2)*b[t]+WLCxi/.{xi->xi[t],eta->eta[t]};
bdot = -2*(xi^2+eta^2)*a[t]+WLCeta/.{xi->xi[t],eta->eta[t]};
timedot = xi[t]^2+eta[t]^2; (* unregularized time variable could be useful *)
LCRegularizedvf = {xidot, etadot, adot, bdot};

LCRegularizedode={xi'[t]==xidot,eta'[t]==etadot, a'[t]==adot,b'[t]==bdot,time'[t]==timedot};

MakeLCRegularizedOrbit[{xi0_,eta0_,a0_,b0_},mu0_,t0_,t1_]:=Module[{h0},
h0 = LeviCivitaEnergy[{xi0,eta0,a0,b0},mu0];
NDSolveValue[(LCRegularizedode/.{mu->mu0,h->h0})~Join~{xi[0]==xi0,eta[0]==eta0,a[0]==a0,b[0]==b0,time[0]==0},{xi[t],eta[t],a[t],b[t],time[t]},{t,t0,t1}]]

Since the motivation for this regularization is to study orbit near second primary, it’s reasonable to use Hill’s regions to help choose interesting initial conditions.

(WLeviCivita/.{mu->mu0,h->h0})

MakeLCHillPlot[mu0_,h0_,options___] :=
RegionPlot[(WLeviCivita/.{mu->mu0,h->h0})≥ 0,{xi,-2,2},{eta,-2,2},PlotStyle->{LightBlue},PlotPoints->30,options]

lchillplot=MakeLCHillPlot[1/3,-2.0,PlotRange->{{-1,1},{-1,1}},Axes->True,AxesLabel->{ξ,η}]

We need a function to set the energy in the Levi-Civita variables.

SetLCEnergy[{xi0_,eta0_,a0_,b0_},mu0_,h0_]:=Module[{vnormsq,a1,b1},
vnormsq = 2*WLeviCivita/.{xi->xi0,eta->eta0,mu->mu0,h->h0}; (*the desired velocity norm squared *)
If [vnormsq<0, Print["Position not in Hill's region"];Return[{xi0,eta0,a0,b0}]];
{a1,b1}={a0,b0}*Sqrt[Abs[vnormsq]]/Sqrt[a0^2+b0^2];
{xi0,eta0,a1,b1}]

Here is a retrograde orbit around the primary starting on the positive eta axis and moving right orthogonally

mu0=1/3.;
xi0=0;eta0=0.2;
a0=1;b0=0; (*this will be normalized to the right energy below *)
h0=-2.0;
{xi0,eta0,a0,b0}= SetLCEnergy[{xi0,eta0,a0,b0},mu0,h0]
orb3=MakeLCRegularizedOrbit[{xi0,eta0,a0,b0},mu0,0,100]
Show[lchillplot,ParametricPlot[orb3[[{1,2}]],{t,0,100}]]

Finding the retrograde periodic orbit by “shooting” and an isolating block

For energy levels where there is a disklike Hill’s region around P2, there is a simple, symmetric periodic orbit which moves exactly once around P2 before closing up.  It can be found by a shooting argument.  Consider all possible initial conditions of the form (xi,eta,a,b) = (0,eta,a,0) with eta>0 and a>0 and which lie in the Hill’s region.  The value of a will  be uniquely determined by eta and so this is a one-parameter family of initial condition.  Among these there is one for which the solution crosses the first quadrant and meet the positive xi axis orthogonally.   The demonstrations in this section provide a quick experimental “proof” of this.  By moving the slider one can find the initial point with the right behavior.

First use Mathematica’s WhenEvent option to calculate an orbit segment up to the moment when it either start to move vertically or else hits the xi axis.   This corresponds to the exit criteria  for Conley’s proof based on an isolating block.  Our WhenEvent will be that  either eta = 0 or a=0 and at  that time,  NDSolve will stop numerical integration and return the  solutions and the final time.

whenevent = WhenEvent[{eta[t]==0,a[t]==0} , {Sow[t],"StopIntegration"}];
MakeLCOrbitToBlockExit[eta0_,mu0_,h0_,tmax_]:= Module[{xi1,eta1,a1,b1,sol,t1},
{xi1,eta1,a1,b1}=SetLCEnergy[{0,eta0,1,0},mu0,h0];
{sol,t1}=Reap@NDSolveValue[(LCRegularizedode/.{mu->mu0,h->h0})~Join~{xi[0]==xi1,eta[0]==eta1,a[0]==a1,b[0]==b1,time[0]==0}~Join~{whenevent},{xi[t],eta[t],a[t],b[t],time[t]},{t,0,tmax}];
{sol,t1}]

The next function just plots the orbit segment inside the Hill’s region

lchillplot2=MakeLCHillPlot[1/3,-2.0,PlotRange->{{-.65,.65},{-.65,.65}},Axes->True,AxesLabel->{ξ,η}];MakeLCShootingPlot[eta0_]:=Module[{orb,t1},
{orb, t1}= MakeLCOrbitToBlockExit[eta0,1/3,-2.0,10];
Show[lchillplot2,ParametricPlot[orb[[{1,2}]],{t,0,t1[[1,1]]}]]]

By moving the slider, you can watch how the orbit segments inside Conley’s isolating block change as the initial point is changed.  If the initial eta is large (near the zero velocity curve) the solution become vertical before it reaches the xi axis.  If the initial eta is near 0 the solution  reaches the xi axis but before it becomes vertical.  The desired orbit is the one  which just reaches the xi axis at the same time it becomes vertical, somewhere around eta = 0.36 in this case.

Manipulate[MakeLCShootingPlot[eta0],{{eta0,0.3},0,0.6}]

It’s possible to find the required initial condition for the periodic  orbit  more precisely.  This time we want to follow orbits all the way the xi axis by changing the WhenEvent (orbit close to  the  right  one  will all eventually reach the axis).  Then  we want to read off the final value of  a = axidot.  This will be zero for the correct orbit.  Here is a function which reads off the final value of a.

whenevent2 = WhenEvent[{eta[t]==0} , {Sow[a[t]],"StopIntegration"}];
FindFinalaValue[eta0_,mu0_,h0_,tmax_]:= Module[{xi1,eta1,a1,b1,sol,t1},
{xi1,eta1,a1,b1}=SetLCEnergy[{0,eta0,1,0},mu0,h0];
a1=Reap@NDSolveValue[(LCRegularizedode/.{mu->mu0,h->h0})~Join~{xi[0]==xi1,eta[0]==eta1,a[0]==a1,b[0]==b1,time[0]==0}~Join~{whenevent2},{xi[t],eta[t],a[t],b[t],time[t]},{t,0,tmax}];
a1[[2,1,1]]]

This plots the final a value as a function of initial eta.  Clearly there is a zero between 0.35 and 0.36.

FindFinalaValue[0.37,1/3,-2,10]
Plot[FindFinalaValue[eta0,1/3,-2,10],{eta0,0.35,0.37}]

We can use bisection to home in on the exact value.  The next function bisects and interval where a function has opposite signs  and returns another such interval, half the size.

Bisect[f_,{x0_,x1_}]:= Module[{x2},
x2=(x0+x1)/2;
If[f[x0]*f[x2]<=0,{x0,x2},If[f[x2]*f[x1]<=0,{x2,x1},Print["Endpoints have the same sign"];{x0,x1}]]]

This wil be applied repeatedly to the function  FindFinalaValue[ #,1/2,-2,10]&.

Bisect[FindFinalaValue[#,1/3,-2,10]&,{0.35,0.36}]
{eta0,eta1}=Nest[Bisect[FindFinalaValue[#,1/3,-2,10]&,#]&,{0.35,0.36},20]

Now use eta0 as an initial value for the periodic orbit. The period is 4 *t1 where t1 is the time to reach the xi axis.

{orbseg,t1} = MakeLCOrbitToBlockExit[eta0,1/3,-2.0,10]
{xi1,eta1,a1,b1} = SetLCEnergy[{0,eta0,1,0},1/3,-2]
periodicorb = MakeLCRegularizedOrbit[{xi1,eta1,a1,b1},1/3,0,4*t1[[1]]]

pp1= ParametricPlot[periodicorb[[{1,2}]],{t,0,t1[[1,1]]},PlotStyle->{Thicker,Red}];
pp2= ParametricPlot[periodicorb[[{1,2}]],{t,0,4*t1[[1,1]]},PlotStyle->Blue];
Show[lchillplot2,pp2,pp1,Graphics[{Green,PointSize[0.03],Point[{0,0}]}]]