Chapter 3, Section 3.8

Here we start by making use of the function Stability (see Sect.3.7.nb) which returns, among other things, an analytical solution for motion around the chosen equilibrium point.

Stability[μ_, ptid_Integer, {X0_, Y0_}] :=
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];
Print["Mass ratio = ", μ,"\n"];
Print["Lagrangian Point ", ptid,"\n"];
Print["{x0, y0} = ", {x,y},"\n"];
If[ptid == 1 || ptid == 2 || ptid ==3,
   abar = (1-μ)/r1^3 + μ/r2^3;
   ceqq = λ^4 + (2-abar)*λ^2 + (1+abar-2*abar^2);
   soln = Chop[λ /. NSolve[ceqq == 0, {λ}]]];
If[ptid == 4 || ptid == 5,
   ceqq = λ^4 + λ^2 + (27/4)*μ2*(1-μ2);
   soln = Chop[λ /. NSolve[ceqq == 0, {λ}]]];
abar = μ1/r1^3 + μ2/r2^3;
cbar = 3*(μ1*(x+μ2)/r1^5 + μ2*(x-μ1)/r2^5)*y;
dbar = 3*(μ1*(x+μ2)^2/r1^5 + μ2*(x-μ1)^2/r2^5);
uxx = 1 - abar + dbar;
uxy = cbar;
Print["Eigenvalues :"];
Print[ColumnForm[soln],"\n"];
Do[s[i] = soln[[i]], {i,1,4}];
Do[b[i] = (s[i]^2 - uxx)*a[i]/(2*s[i] + uxy), {i,1,4}];
eq1 = a[1] + a[2] + a[3] + a[4] == X0;
eq2 = b[1] + b[2] + b[3] + b[4] == Y0;
eq3 = s[1]*a[1] + s[2]*a[2] + s[3]*a[3] + s[4]*a[4] == 0;
eq4 = s[1]*b[1] + s[2]*b[2] + s[3]*b[3] + s[4]*b[4] == 0;
sol = Flatten[Chop[NSolve[
         {eq1,eq2,eq3,eq4},{a[1], a[2], a[3], a[4]}]]];
asol = {a[1], a[2], a[3], a[4]} /. sol;
bsol = {b[1], b[2], b[3], b[4]} /. sol;
Print["\nValues of α bar: "];
Print[ColumnForm[asol]];
Print["\nValues of β bar: "];
Print[ColumnForm[bsol]];
xx[t_] := 
Chop[ComplexExpand[
 ({a[1], a[2], a[3], a[4]} . 
    {Exp[s[1]*t], Exp[s[2]*t], Exp[s[3]*t], Exp[s[4]*t]}) /.
sol]];
yy[t_] := 
Chop[ComplexExpand[
 ({b[1], b[2], b[3], b[4]} . 
    {Exp[s[1]*t], Exp[s[2]*t], Exp[s[3]*t], Exp[s[4]*t]}) /.
sol]];
Print["\nX(t) = "];
Print[xx[t]];
Print["\nY(t) = "];
Print[yy[t],"\n"]
]

We are going to study the solution around the L4 point for a mass ratio of 0.01.

Stability[0.01, 4, {0.00001,0.00001}]

The X(t) and Y(t) solutions given above were stored internally in Stability as xx[t] and yy[t].  We can illustrate the fact that there are two fundamental motions by carrying out a rotation of 30° clockwise about the equilibrium point.  First we define the rotation matrix, rot. (see Eq.(3.152).)

rot = N[{{Cos[30°], Sin[-30°]},
         {Sin[30°], Cos[30°]}}];

MatrixForm[rot]

Then we apply rot to the functions xx[t] and yy[t], "chopping" to remove any small imaginary parts.  (cf. Eq.(3.153).)

ctranxy = Chop[Expand[rot . {xx[t], yy[t]}],0.00001]

Note that the ratio of the magnutude of the two short period terms is approximately 2, implying that the motion around the epicycle is approximately that of a 2:1 centred ellipse.

0.00009853667909211485/0.00004861047011093332

Plot the trajectory that results from the analytical solution for the motion around the equilibrium point.

traj=
ParametricPlot[{xx[t],yy[t]},{t,0,25Pi},
PlotPoints->50,
PlotRange->{{-0.0005,0.0005},{-0.00025,0.00025}},
AspectRatio->0.5,Axes->None]

Now define a function czvn that returns the value of the Jacobi constant of the zero velocity curve at the point (x,y) where x and y are the coordinates referred to the L4 equilibrium point.  (See Eq.(3.39).)

czvn[x_,y_,μ_] :=
Module[{xn = x + 1/2 - μ, yn = y + Sqrt[3]/2}, 
xn^2+yn^2+
2*((1-μ)/Sqrt[(xn+μ)^2+yn^2] + μ/Sqrt[(xn-1+μ)^2+yn^2])
];

Now we define a function jacobin that gives the value of the Jacobi constant for a given position and velocity in the rotating frame and for a given mass ratio, μ.  Note that the position (x,y) and velocity (xd,yd) are expressed with respect to the L4 point.

jacobin[{x_, y_, xd_, yd_}, μ_] := 
Module[{xb = x+1/2-μ, yb = y + N[Sqrt[3]/2]},
xb^2+yb^2-xd^2-yd^2+ 
2*(μ/Sqrt[(-1+μ+xb)^2+yb^2]+(1-μ)/Sqrt[(μ+xb)^2+yb^2])];

Calculate the value of the Jacobi constant for our trajectory by evaluating it at the starting point using jacobin.

jcon = jacobin[{xx[0], yy[0], 0, 0}, 0.01]

Plot the zero velocity curve for this value of the Jacobi constant.

zvcurve=
ContourPlot[czvn[x,y,0.01],
{x,-0.0004,0.0004},{y,-0.0002,0.0002},PlotPoints->100,
AspectRatio->Automatic,Axes->None,
Contours->{jcon},ContourShading->False];

Plot the zero velocity curve and the trajectory on the same plot.  Note how the trajectory can touch but never cross the zero velocity curve.

trajcurve=
Show[zvcurve,traj,Graphics[Point[{0,0}]],
AspectRatio->Automatic]

Remind ourselves of the solution for the compound motion.  Note the two frequencies in the motion.

ctranxy

Let us choose a particular time at which to illustrate the solution.

tt = 0;

We can calculate the position of the guiding centre at this time by isolating just the long period part of the motion.

centre1 = 
{0.00035373*Sin[0.268348*t], 0.0000622707*Cos[0.268348*t]} /. t -> tt

Now we reverse the effects of the previous rotation by creating a rotation of 30° in the opposite direction.

rot = N[{{Cos[30°], Sin[30°]},
         {Sin[-30°], Cos[30°]}}]

Construct a dashed-line plot showing the epicyclic motion at this time by adding centre1 to the short period terms in the solution (and reversing the previous rotation).

ellipse1 = 
ParametricPlot[
 rot.(centre1 + 
  {-0.0000985367*Sin[0.963322*t],
   -0.0000486105*Cos[0.963322*t]}),
 {t,0,2Pi/0.963322},
 AspectRatio -> Automatic,
 PlotStyle->Dashing[{0.01,0.01}]]

Do the same to plot the motion of the epicentre by isolating just the long period terms.

ellipse2 = 
ParametricPlot[
 rot.{0.00035373*Sin[0.268348*t], 0.0000622707*Cos[0.268348*t]},
 {t,0,2Pi/0.268348},
 AspectRatio -> Automatic,
 PlotStyle->Dashing[{0.01,0.01}]]

Show both plots together without axes, adding a point to denote the location of the equilibrium position.

ellipses = 
Show[ellipse1,ellipse2,Graphics[Point[{0,0}]],
     PlotRange->{{-0.0005,0.0005},{-0.00025,0.00025}},
     AspectRatio->Automatic,Axes->None]

Finally we put all the plots together to illustrate the type of motion.

Show[traj,ellipses]