Chapter 8, Section 8.12.1, Figure 8.23

The definitions of the following functions are taken from the PhD thesis of Malhotra (1988).

First we define the function ψ.

ψ[δ_] := N[ArcTan[Sqrt[(δ+1)^3-1]]];

Next we define the function Φ3.

Φ3[δ_] := N[(1/2)*(δ+1)*(Cos[ψ[δ]/3]+Sqrt[3]*Sin[ψ[δ]/3])^2];

The functions Φmin and Φmax determine the minimum and maximum values of Φ.  These are given in Eqs.(8.211) and (8.212).

Φmin[δ_] := N[3*(δ+1) - Φ3[δ] - 2*(2*Φ3[δ])^(1/4)];

Φmax[δ_] := N[3*(δ+1) - Φ3[δ] + 2*(2*Φ3[δ])^(1/4)];

The function γ is defined in Eq.(8.215).

γ[δ_] := 
N[ArcSin[(Φmax[δ]+Φmin[δ]-2*Φ3[δ])/(Φmax[δ]-Φmin[δ])]];

The functions Areaint and Areaext are the same as those defined in Eq.(8.213) and (8.214).

Areaint[δ_] := N[Abs[
3*((1/2)*(Φmax[δ] + Φmin[δ] + 2*Φ3[δ] - 2*(δ+1))*
(Pi/2-γ[δ]) - Sqrt[(Φmax[δ]-Φ3[δ])(Φ3[δ]-Φmin[δ])])]];

Areaext[δ_] := N[Abs[
-3*((1/2)*(Φmax[δ] + Φmin[δ] + 2*Φ3[δ] - 2*(δ+1))*
 (Pi/2+γ[δ]) + Sqrt[(Φmax[δ]-Φ3[δ])(Φ3[δ]-Φmin[δ])])]];

The x values of the equilibrium points are the solutions of the cubic equation defined by xval.  Note that this is just Eq.(8.135) with the transformation given in Eq.(8.206).

xval[δ_]:=NSolve[x^3-3(δ+1)x-2 == 0, x];

The function Φval takes the real solutions given by xval and calculates the corresponding values of Φ.

Φval[δ_] := Sqrt[2*Max[Select[x /. xval[δ],(Im[#]==0)&]]]

Now define the Hamiltonian, h1, of the system.

h1[x_,y_,δ_]:=
   -(3/2)*(δ+1)*(x^2+y^2)+(x^2+y^2)^2/4-2*x;

The function xpts returns the values of x for y = 0 that lie on a particular curve of constant Hamiltonian.

xpts[hval_, δ_] := 
Select[x /. NSolve[h1[x,0,δ] == hval, x], (Im[#]==0)&]

The function ymax returns the maximum value of y for particular values of x, the Hamiltonian and δ.

ymax[x_, h_, δ_] :=
Module[{a,b,c},
a = (1/4);
b = (1/2)*x^2 - (3/2)*(δ+1);
c = (1/4)x^4 - (3/2)*(δ+1)*x^2 - 2*x - h;
Sqrt[(-b + Sqrt[b^2-4*a*c])/2/a]]

The function area calculates the area enclosed by a curve defined by the Hamiltonian value h for a particular δ.

area[h_,δ_] := 
Module[{xmin,xmax},
{xmin,xmax} = xpts[h,δ];
2*NIntegrate[y,{x,xmin,xmax},{y,0,ymax[x,h,δ]}]]

The function ΦValue returns the value of Φ for a particular δ and Hamiltonian h.

ΦValue[h_,δ_] := N[area[h,δ]/(2*Pi)]

The function h1Φ returns the Hamiltonian value taking input in terms of Φ, ϕ, and δ.

h1Φ[Φ_,ϕ_,δ_]:=N[-3(δ+1)*Φ + Φ^2 - 2*Sqrt[2*Φ]*Cos[ϕ]];

In the series of plots we are assuming that Φ, the action or alternatively the area enclosed by the trajectory, has a value of 4.  We start with δ = –6 and find a value of the Hamiltonian that gives Φ = 4.

test = ΦValue[65.5905,-6]

Now calculate the value of xmin and xmax for this value of the Hamiltonian.

{xmin, xmax} = xpts[65.5905,-6]

Now calculate the value of Φleft and Φright for these values.

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

Verify that this is indeed the correct value of the Hamiltonian.

h = N[h1Φ[Φright,0,-6]]

Plot the curve associated with this value of the Hamiltonian.  See Fig.8.23a.

p21 =
ContourPlot[h1[x,y,-6],{x,-5,5},{y,-5,5},
   ContourShading->False, AxesOrigin->{0,0},
   Contours->{h},PlotPoints->100]

Now go through the same procedure for δ = 0.

test = Chop[ΦValue[1.99648,0]]

{xmin, xmax} = xpts[1.99648,0]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,0]]

Here we also calculate the value of the Hamiltonian associated with the critical curve.

hcrit = h1Φ[Φmax[0],0,0]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.22b.

p23 =
ContourPlot[h1[x,y,0],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->100]

Now go through the same procedure for δ = 0.485584.

N[Areaext[0.485584]/2/Pi]

h = h1Φ[Φmax[0.485584],0,0.485584]

hcrit = h1Φ[Φmax[0.485584],0,0.485584]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.23c,d,f.

p24 =
ContourPlot[h1[x,y,0.485584],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{hcrit},PlotPoints->150]

Now we redefine ΦValue to look at the capture case.

newΦValue[h_, δ_] :=
Module[{Φ,Φleft,Φright,area,xmin,xmax,temp},
temp = -3*(δ+1)*Φ+Φ^2-h;
{xmin,xmax} = xpts[h,δ];
Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2];
Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2];
Abs[Chop[2*NIntegrate[
(-3*(δ+1)*Φ+2*Φ^2-(1/2)*temp)/Sqrt[8*Φ-temp^2],
{Φ,Φleft,Φright}, SingularityDepth->8,
MaxRecursion->12,WorkingPrecision->20]]]
]

We are now using δ = 4.

N[newΦValue[-51.85,4]/2/Pi]

newΦ = 
N[Areaext[0.485584]/2/Pi]-N[Areaint[0.485584]/2/Pi]

N[Areaext[4]/2/Pi]-N[Areaint[4]/2/Pi]

N[Areaint[4]/2/Pi]

{xmin, xmax} = xpts[-51.85,4]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,4]]

hcrit = h1Φ[Φmax[4],0,4]

Plot the curves associated with these values of the Hamiltonian.  See Fig.8.23e.

p25 =
ContourPlot[h1[x,y,4],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->150]

Now we redefine ΦValue to look at the passage case.

newΦValue[h_, δ_] :=
Module[{Φ,Φleft,Φright,area,xmin,xmax,temp},
temp = -3*(δ+1)*Φ+Φ^2-h;
{xmin,xmax} = Take[xpts[h,δ],{2,3}];
Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2];
Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2];
Abs[Chop[2*NIntegrate[
(-3*(δ+1)*Φ+2*Φ^2-(1/2)*temp)/Sqrt[8*Φ-temp^2],
{Φ,Φleft,Φright}, SingularityDepth->8,
MaxRecursion->12,WorkingPrecision->20]]]
]

Again we are using δ = 4.

N[newΦValue[-11.1,4]/2/Pi]

xpts[-11.1,4]

newΦ = N[Areaint[0.485584]/2/Pi]

N[Areaext[4]/2/Pi]-N[Areaint[4]/2/Pi]

N[Areaint[4]/2/Pi]

{xmin, xmax} = Take[xpts[-11.1,4],{2,3}]

Φleft = If[xmin < 0, -Abs[xmin]^2/2, Abs[xmin]^2/2]

Φright = If[xmax < 0, -Abs[xmax]^2/2, Abs[xmax]^2/2]

h = N[h1Φ[Φright,0,4]]

hcrit = h1Φ[Φmax[4],0,4]

Plot the curves associated with these values of the Hamiltonian.  Note that we are not interested in the outer branch of the non-critical curve.  See Fig.8.23g.

p35 =
ContourPlot[h1[x,y,4],{x,-5,5},{y,-5,5},
   ContourShading->False, 
   Contours->{h,hcrit},PlotPoints->150]