ClearAll[e]
phi[x_,y_]=(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_]=-D[phi[x,y],x]
fy[x_,y_]=-D[phi[x,y],y]
rot[t_]={{Cos[2 Pi t], -Sin[2 Pi t]},{Sin[2 Pi t],Cos[2 Pi t]}}

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

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


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

Simplify[phi[xp,yp,e]]

if moon:
 e=.0121

(* set your guess for initial speed, angle (0 radians is right at Moon), initial x value (initial y is set to zero) *)
v0=XX; theta=XX; x0=XX

solution=NDSolve[{x''[t]==fx[x[t],y[t]]+(4 Pi) y'[t],
y''[t]==fy[x[t],y[t]]-(4 Pi) x'[t],
x[0]==x0,y[0]==0,
x'[0]==v0 Cos[theta],y'[0]==v0 Sin[theta]},
{x,y},{t,0,1},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]+e,y[t]} /. solution],{t,0,.3},AspectRatio->Automatic]


(* if doing Lenz stuff *)
vx[t_]=x'[t]-y[t]2 Pi /. First[solution]
vy[t_]=y'[t]+x[t]2 Pi /. First[solution]
xe[t_]=x[t]+e /. 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]
