phi[x_,y_,e_]=(2 Pi)^2 ( -(1-e)/Sqrt[(x+e)^2+y^2] -e/Sqrt[(x-(1-e))^2+y^2] -
     (1/2)(x^2+y^2)    )
fx[x_,y_,e_]=-D[phi[x,y,e],x]
fy[x_,y_,e_]=-D[phi[x,y,e],y]

Plot3D[phi[x,y,.1],{x,-1.5,1.5},{y,-1.5,1.5},
PlotPoints->50,PlotRange->{-80,-50},ViewPoint->{2,-2,.6}]

ContourPlot[phi[x,y,.1],{x,-1.5,1.5},{y,-1.2,1.2},
Contours->10,PlotPoints->50,PlotRange->{-80,-50},AspectRatio->Automatic]


xp=-e+Cos[Pi/3]
yp=Sin[Pi/3]
Simplify[fx[xp,yp,e]]
Simplify[fy[xp,yp,e]]

stability of L4: e=.01 & e=.04
solution=NDSolve[{x''[t]==fx[x[t],y[t],.01]+(4 Pi) y'[t], 
    y''[t]==fy[x[t],y[t],.01]-(4 Pi) x'[t], 
    x[0]==-.01+Cos[Pi/3]+.01, y[0]==Sin[Pi/3], x'[0]==0, y'[0]==0}, 
    {x,y},{t,0,10}]
ParametricPlot[Evaluate[{x[t],y[t]} /. solution],{t,0,10}]
ContourPlot[phi[x,y,.01],{x,0,1.3},{y,0,1.3},ContourShading -> False,
Contours->20,PlotPoints->50,PlotRange->{-65,-59},AspectRatio->Automatic]
Show[%,%%,AspectRatio->Automatic]

solution=NDSolve[{x''[t]==fx[x[t],y[t],.04]+(4 Pi) y'[t], y''[t]==fy[x[t],y[t],.04]-(4 Pi) x'[t], 
    x[0]==-.04+Cos[Pi/3]+.001, y[0]==Sin[Pi/3],
    x'[0]==0,         y'[0]==0}, {x,y},{t,0,10}]
ParametricPlot[Evaluate[{x[t],y[t]} /. solution],{t,0,10}]
ContourPlot[phi[x,y,.01],{x,0,1.3},{y,0,1.3},ContourShading -> False,
Contours->20,PlotPoints->50,PlotRange->{-65,-59},AspectRatio->Automatic]
Show[%,%%,AspectRatio->Automatic]


moon shot:
escape velocity:
Sqrt[2*( phi[xp,yp,e] -phi[.1,0,e])] /. e->.0121
Out[19]= 24.0692

v0=25.; theta=.42

solution=NDSolve[{x''[t]==fx[x[t],y[t],.0121]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t],.0121]-(4 Pi) x'[t],
x[0]==.1,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,.6},MaxSteps-> 50000,AccuracyGoal -> 12]

ParametricPlot[Evaluate[{x[t],y[t]}/. solution],{t,0,.4},
AspectRatio->Automatic,PlotRange->All,PlotPoints->200]

rot[t_]={{Cos[2 Pi t], -Sin[2 Pi t]},{Sin[2 Pi t],Cos[2 Pi t]}}

ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,0,.4},
AspectRatio->Automatic]

orbits:
v0=.6*25./Sqrt[2]; theta=Pi/2

solution=NDSolve[{x''[t]==fx[x[t],y[t],.0121]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t],.0121]-(4 Pi) x'[t],
x[0]==.1,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,.6},MaxSteps-> 50000,AccuracyGoal -> 12]

ParametricPlot[Evaluate[{x[t],y[t]}/. solution],{t,0,.4},
AspectRatio->Automatic,PlotRange->All,PlotPoints->200]

ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,0,.4},
AspectRatio->Automatic]

solution=NDSolve[{x''[t]==fx[x[t],y[t],.0121]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t],.0121]-(4 Pi) x'[t],
x[0]==.1,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,4},MaxSteps-> 50000,AccuracyGoal -> 16]

ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,3.81,3.8283},
AspectRatio->Automatic]

energy conservation:
Plot[Evaluate[(x'[t]^2+y'[t]^2)/2+phi[x[t],y[t],.0121]  /. solution],{t,0,3.8}]

lenz:
vx[t_]=x'[t]-y[t]2 Pi /. First[solution]
vy[t_]=y'[t]+x[t]2 Pi /. First[solution]
xe[t_]=x[t]+.0121 /. First[solution]

lenz=rot[t].({vy[t](y[t]vx[t]-xe[t]vy[t]),vx[t](xe[t]vy[t]-y[t]vx[t])}/(2 Pi)^2 + 
{xe[t],y[t]}/Sqrt[xe[t]^2+y[t]^2]) /. First[solution]
ParametricPlot[Evaluate[lenz],{t,0,3}]
ParametricPlot[Evaluate[lenz],{t,0,.02}]
ParametricPlot[Evaluate[lenz],{t,2.98,3}]

slingshot:
v0=25.; theta=.405

solution=NDSolve[{x''[t]==fx[x[t],y[t],.0121]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t],.0121]-(4 Pi) x'[t],
x[0]==.1,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,.6},MaxSteps-> 50000,AccuracyGoal -> 12]

ParametricPlot[Evaluate[{x[t],y[t]}/. solution],{t,0,.4},
AspectRatio->Automatic,PlotRange->All,PlotPoints->200]
ParametricPlot[Evaluate[{x[t],y[t]}/. solution],{t,.099,.102}, 
AspectRatio->Automatic,PlotRange->All,PlotPoints->200]         

ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,0,.4},
AspectRatio->Automatic]
ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,.099,.102},
AspectRatio->Automatic]

fixed-frame energy gain:
vx[t_]=x'[t]-y[t]2 Pi /. First[solution]
vy[t_]=y'[t]+x[t]2 Pi /. First[solution]

Plot[Evaluate[(vx[t]^2+vy[t]^2)/2+phi[x[t],y[t],.0121] +(2 Pi)^2 (1/2)(x[t]^2+y[t]^2) /. solution],{t,0,.4}]

compare:moon shot

v0=25.; theta=.42

solution=NDSolve[{x''[t]==fx[x[t],y[t],.0121]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t],.0121]-(4 Pi) x'[t],
x[0]==.1,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,.6},MaxSteps-> 50000,AccuracyGoal -> 12]
ParametricPlot[Evaluate[rot[t].{x[t]+.0121,y[t]} /. solution],{t,0,.4},
AspectRatio->Automatic]

vx[t_]=x'[t]-y[t]2 Pi /. First[solution]
vy[t_]=y'[t]+x[t]2 Pi /. First[solution]

fixed-frame energy loss:
Plot[Evaluate[(vx[t]^2+vy[t]^2)/2+phi[x[t],y[t],.0121] +(2 Pi)^2 (1/2)(x[t]^2+y[t]^2) /. solution],{t,0,.2}]

