Central Force Problem

This notebook has some programs for finding orbits of central force problems and plotting them in interesting ways.

Phase Portraits for the Reduced  Central Force Problem

A central force problem in the plane is determined by giving a potential function F(r) (the potential energy is -F(r)).  Here are some examples.  For simplicity, the mass parameter is set to m=1 for the power law.

FPowerLaw[r_,a_]:= 1/r^a  (* Kepler is a=1.  Mass is m=1 *)
FOblatePlanet[r_,delta_,m_]:=m/r + delta*m/(2*r^3)

Fixing an angular momentum c gives a reduced system in the (r,w) plane where  w = rdot.  The phase curves are level curves of the energy function.  The amended potential energies are also defined here.

EnergyPowerLaw[r_,w_,c_,a_]:=    w^2/2 +c^2/(2*r^2) -FPowerLaw[r,a]  
EnergyOblatePlanet[r_,w_,c_,delta_,m_]:=    w^2/2 +c^2/(2*r^2) -FOblatePlanet[r,delta,m]  
VPowerLaw[r_,c_,a_]:= c^2/(2*r^2) -FPowerLaw[r,a] 
VOblatePlanet[r_,c_,delta_,m_]:= c^2/(2*r^2) -FOblatePlanet[r,delta,m]

Here is the amended potential energy and the energy level plot for the Kepler problem and the power law 1/r^3 for angular momentum c=Sqrt[2].

p1= Plot[VPowerLaw[r,Sqrt[2],1],{r,0,10}, PlotPoints->40,AxesLabel->{r,"V(r)"},PlotLabel->"α=1"];
en1=ContourPlot[EnergyPowerLaw[r,w,Sqrt[2],1],{r,0,6},{w,-1,1},ContourShading->False,Contours->20, PlotPoints->30];

p3= Plot[VPowerLaw[r,Sqrt[2],3],{r,0,10}, PlotPoints->40,AxesLabel->{r,"V(r)"},PlotLabel->"α=3",PlotRange->{-0.1,0.25}];
en3=ContourPlot[EnergyPowerLaw[r,w,Sqrt[2],3],{r,0,6},{w,-1,1},ContourShading->False,Contours->20, PlotPoints->30];

Show[GraphicsGrid[{{p1,en1}}]]
Show[GraphicsGrid[{{p3,en3}}]]

The power law potentials for a < 2 look like the Kepler one and for a >2 they resemble the a = 3 case.

Manipulate[Plot[VPowerLaw[r,Sqrt[2],a],{r,0,10}, PlotPoints->40,PlotLabel->"V(r,α)"],{a,0,4}]

Solving the Differential Equations

First set up the differential equations in Mathematica notation.  Here Fr is the derivative of the potential.

FrPowerLaw[r_,a_] = D[FPowerLaw[r,a],r]
FrOblatePlanet[r_,delta_,m_] = D[FOblatePlanet[r,delta,m],r]

ODEPowerLaw[a_] = With[{r=Sqrt[x[t]^2+y[t]^2]},{x'[t]== u[t],y'[t]==v[t],u'[t]==x[t]*FrPowerLaw[r,a]/r,v'[t]==y[t]*FrPowerLaw[r,a]/r}]

ODEOblatePlanet[delta_,m_] = With[{r=Sqrt[x[t]^2+y[t]^2]},{x'[t]== u[t],y'[t]==v[t],u'[t]==x[t]*FrOblatePlanet[r,delta,m]/r,v'[t]==y[t]*FrOblatePlanet[r,delta,m]/r}]//Simplify

Next we need to choose some initial conditions.  For the Kepler problem, it’s reasonable to do this using orbital elements a, e, varpi.   Another way, which works for other problems is  just giving the initial {x0,y0,u0,v0}.  Alternatively, you can specify a reduced initial condition {r0,w0} and an angular momentum, c.

KeplerIC[a_,ecc_,varpi_]:= Module[{x0,y0,u0,v0},
{x0,y0}=a*(1-ecc)*{Cos[varpi],Sin[varpi]};
{u0,v0}=Sqrt[(1+ecc)/Abs[a*(1-ecc)]]*{-Sin[varpi],Cos[varpi]};
{x[0]==x0,y[0]==y0,u[0]==u0,v[0]==v0}
]

MakeIC[{x0_,y0_,u0_,v0_}]:= {x[0]==x0,y[0]==y0,u[0]==u0,v[0]==v0}

MakeIC[{r0_,w0_},c_]:= {x[0]==r0,y[0]==0,u[0]==w0,v[0]==c/r0}

Now use NDSolve to compute orbits.

MakeKeplerOrbit[a_,ecc_,varpi_,tmax_]:=NDSolveValue[ODEPowerLaw[1]~Join~KeplerIC[a,ecc,varpi],{x[t],y[t],u[t],v[t]},{t,-tmax,tmax}]

MakeKeplerOrbit[{x0_,y0_,u0_,v0_},tmax_]:=NDSolveValue[ODEPowerLaw[1]~Join~MakeIC[{x0,y0,u0,v0}],{x[t],y[t],u[t],v[t]},{t,-tmax,tmax}]

MakePowerLawOrbit[{r0_,w0_},c_,alpha_,tmax_]:=NDSolveValue[ODEPowerLaw[alpha]~Join~MakeIC[{r0,w0},c],{x[t],y[t],u[t],v[t]},{t,-tmax,tmax}]

Kepler Orbits

orb= MakeKeplerOrbit[1,0.7,Pi/4,10]
ParametricPlot[orb[[{1,2}]],{t,-10,10},PlotLabel->"Orbit"];
ParametricPlot[orb[[{3,4}]],{t,-10,10},PlotLabel->"Hodograph",PlotStyle->Green];
GraphicsGrid[{{%%,%}}]

Here is an interactive plotter.  First make a function which draws an orbit an puts a dot on it corresponding to position at time t0.

MakeOrbitPlot[orb_,t0_,tmax_,plotopts___]:= Module[{pp,q,v},
pp= ParametricPlot[orb[[{1,2}]], {t,-tmax,tmax},PlotPoints->Floor[100*tmax],plotopts];
q= Graphics[{Red,PointSize[0.03],Point[orb[[{1,2}]]/.t->t0]}];
v= Graphics[{Black,Arrow[{orb[[{1,2}]],orb[[{1,2}]]+0.3*orb[[{3,4}]]}/.t->t0]}];
Show[pp,q]]

MakeOrbitPlotWithExtra[orb_,t0_,tmax_,extra_,plotopts___]:= Module[{pp,q,v},
pp= ParametricPlot[orb[[{1,2}]], {t,-tmax,tmax},PlotPoints->Floor[100*tmax],plotopts];
q= Graphics[{Red,PointSize[0.03],Point[orb[[{1,2}]]/.t->t0]}];
v= Graphics[{Black,Arrow[{orb[[{1,2}]],orb[[{1,2}]]+0.3*orb[[{3,4}]]}/.t->t0]}];
Show[pp,q,extra]]

tmax0 = 10;
Manipulate[MakeOrbitPlot[MakeKeplerOrbit[a,ecc,varpi,tmax0],t0,tmax0,PlotRange->2*{{-a,a},{-a,a}}],{{a,1},0,2},{{ecc,0.5},0,2},{{varpi,0} ,-Pi,Pi},{{t0,0},-tmax0,tmax0}]

Here is another interactive plotter where you can select the initial position and velocity directly.

Manipulate[MakeOrbitPlotWithExtra[MakeKeplerOrbit[q~Join~v,tmax0],t0,tmax0,Graphics[Arrow[{q,q+v}]],PlotRange->{{-size,size},{-size,size}}],{{t0,0},-tmax0,tmax0,ImageSize->Large},{{tmax0,10},10,100,ControlPlacement->Right,ImageSize->Small},
{{size,3},0.1,10,ControlPlacement->Right,ImageSize->Small},{{q,{1,0}},Locator},{{q,{1,0}},InputField,ControlPlacement->Right},{{v,{0,1}},InputField,ControlPlacement->Right},{{v,{0,1}},{-1,-1},{1,1},Appearance->"Labeled",ControlPlacement->Right,ImageSize->{150,150}}]

Power Law Orbits

Now for some power law orbits.  To help choose good initial condition, fix and angular momentum and look at the reduced phase portrait.

alpha = 1.5;
c = 1;
ContourPlot[EnergyPowerLaw[r,w,c,alpha],{r,0,3},{w,-1,1},ContourShading->False,Contours->20, PlotPoints->30]

r0=1;
w0=0;
tmax1=20;
size1=1.5;
orb1=MakePowerLawOrbit[{r0,w0},c,alpha,tmax1]
Manipulate[MakeOrbitPlot[orb1,t0,tmax1,PlotRange->{{-size1,size1},{-size1,size1}}],{{t0,0},-tmax1,tmax1}]

Invariant Tori and Rotation Number

A 3D plot of the power law orbits shows the invariant tori.   Since there are really 4 dimensions, you have to defined some third coordinate cooked up from q and v.

r0=1;
w0=0;
tmax0=100;
size=1.5;
orb2=MakePowerLawOrbit[{r0,w0},c,alpha,tmax0]
ParametricPlot3D[{orb2[[1]],orb2[[2]],0.3*orb2[[3]]},{t,-tmax0,tmax0},Axes->False, Boxed->False]

ParametricPlot3D[{orb2[[1]],orb2[[2]],orb2[[{1,2}]].orb2[[{3,4}]]/Norm[orb2[[{1,2}]]]},{t,-tmax0,tmax0},Axes->False, PlotPoints->1000,Boxed->False,ColorFunction->Function[{x,y,z},0.7*z]] (* radial velocity *)

To calculate rotation number of the torus map, we will integrate orbits over half of a period of the radial oscillation and see what happens to the angle.  Mathematica provides a way to stop integration when an Event occurs.  The event will be that the radial velocity w is zero.

MakePowerLawOrbitToWZero[{r0_,w0_},c_,alpha_,tmax_]:=NDSolveValue[ODEPowerLaw[alpha]~Join~MakeIC[{r0,w0},c],{x[t],y[t],u[t],v[t]},{t,0,tmax},
Method->{ "EventLocator","Event"->x[t]*u[t]+y[t]*v[t] }]

r0=1;
w0=0;
tmax0=100;
orb=MakePowerLawOrbitToWZero[{r0,w0},c,alpha,tmax0]

To find the change in theta, we need the last point on this orbit segment, so we need to find the final time.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
InterpolatingFunctionDomain[orb[[1,0]]]
%[[1,2]]

Put all this into a function which reads off the final point on the orbit segment and finds the angle.

FinalAngle[{r0_,w0_},c_,alpha_,tmax_]:=Module[{t1,q1,orb},
orb=MakePowerLawOrbitToWZero[{r0,w0},c,alpha,tmax];
t1=InterpolatingFunctionDomain[orb[[1,0]]][[1,2]];
q1= orb[[{1,2}]]/.t->t1;
Mod[ArcTan[q1[[1]],q1[[2]]],2*Pi]]

r0=1;
w0=0;
tmax0=100;
size=1.5;
ang=FinalAngle[{r0,w0},c,alpha,tmax0]/(2*Pi)

This plot shows that the rotation number varies with the energy for fixed c=1 and power law alpha = 1.5.  So there will be rational and irrational values.

Plot[FinalAngle[{r0,0},1,1.5,100]/(2*Pi),{r0,0.5,2}]

Oblate Planet Orbits

Now for some planar orbits of a satellite around an oblate planet, say the Earth.  It’s a central force problem with potential  m/r + delta*m/(2*r^3) where delta>0 is small.  For time in days and Earth radius 1, the mass is 11468.

delta0 = 0.1; (*this is much larger than the real value of Earth of about 0.001 *)
mearth=11468;
Manipulate[Plot[VOblatePlanet[r,c,delta0,mearth],{r,0,10}, PlotPoints->40,AxesLabel->{r,"V(r)"},PlotLabel->δ==delta0],{c,100,1000}]

Choose a nice angular momentum and look for some satellite orbits

c=150;
ContourPlot[EnergyOblatePlanet[r,w,150,delta0,mearth],{r,0,6},{w,-50,50},ContourShading->False,Contours->20, PlotPoints->30]

MakeOblatePlanetOrbit[{r0_,w0_},c_,delta_,m_,tmax_]:=NDSolveValue[ODEOblatePlanet[delta,m]~Join~MakeIC[{r0,w0},c],{x[t],y[t],u[t],v[t]},{t,-tmax,tmax}]

The orbit shows a prograde precession.

r0=3.0;
w0=0;
tmax2=5;
size2=5;
orb2=MakeOblatePlanetOrbit[{r0,w0},150,delta0,mearth,tmax2]
Manipulate[MakeOrbitPlot[orb2,t0,tmax2,PlotRange->{{-size2,size2},{-size2,size2}}],{{t0,0},-tmax2,tmax2}]

c0=250;
ContourPlot[EnergyOblatePlanet[r,w,c0,delta0,mearth],{r,0,12},{w,-50,50},ContourShading->False,Contours->40, PlotPoints->30]

r0=7;
w0=0;
tmax3=3;
size3=10;
orb3=MakeOblatePlanetOrbit[{r0,w0},c0,delta0,mearth,tmax3]
Manipulate[MakeOrbitPlot[orb3,t0,tmax3,PlotRange->{{-size3,size3},{-size3,size3}}],{{t0,0},-tmax3,tmax3}]