Restricted TBP - evolution

Parameters and  ODEs

Clear["Global`*"]
m=0.0123;
a=(m/3)^(1/3)   (* Hill radius *)

r1=(Sqrt[(x[t]+m)^2+y[t]^2]); r2=(Sqrt[(x[t]-1+m)^2+y[t]^2]);
deq1=x''[t]-2 y'[t]-x[t]==-(1-m)(x[t]+m)/r1^3-m*(x[t]-1+m)/r2^3;
deq2=y''[t]+2x'[t]-y[t]==-y[t]((1-m)/r1^3+m/r2^3);
ivp={deq1,deq2,x[0]==x0, y[0]==y0,x'[0]==vx0, y'[0]==vy0}

Lagrange Points

L4, 5  (+, -)

x04=1/2-m; y04=(Sqrt[3])/2;   
Ci=-0.5*(3-m)
Lpoint4=ListPlot[{{x04,y04}},PlotStyle->{Red,PointSize[Large]}];
Lpoint5=ListPlot[{{x04,-y04}},PlotStyle->{Red,PointSize[Large]}];

L1

x01=1-m-a+a^2/3; y01=0;   
Ci=-0.5*(3+3^(4/3)m^(2/3)-10m/3)
Lpoint1=ListPlot[{{x01,y01}},PlotStyle->{Red,PointSize[Large]}];

L2

x02=1-m+a+a^2/3; y02=0;   
Ci=-0.5*(3+3^(4/3)m^(2/3)-14m/3)
Lpoint2=ListPlot[{{x02,y02}},PlotStyle->{Red,PointSize[Large]}];

L3

x03=-1-m+7m/(12(1-m)); y03=0;   
Ci=-0.5(3+m)
Lpoint3=ListPlot[{{x03,y03}},PlotStyle->{Red,PointSize[Large]}];

Evolution

Initial conditions

x0=x04+0.0065; 
y0=y04+0.0065; 
vx0=0; vy0=0;
tmax=100;

R1=r1/.{x[t]->x0,y[t]->y0};
R2=r2/.{x[t]->x0,y[t]->y0};
Cj=(1/2)(vx0^2+vy0^2) -((1-m)/R1)-(m/R2)-(1/2)(x0^2+y0^2)

Numerical integration

sol=NDSolve[ivp,{x,y},{t,0,tmax}];
xx=x[t]/.sol[[1]]
yy=y[t]/.sol[[1]]

Orbit

orbit=ParametricPlot[{xx,yy},{t,0,tmax}]

Show[{orbit, Lpoint1,Lpoint2,Lpoint3,Lpoint4,Lpoint5}, PlotRange->All, AxesOrigin->{0,0}]

Orbit in Inertial Frame

XX=xx*Cos[t]-yy*Sin[t];
YY=xx*Sin[t]+yy*Cos[t];
ParametricPlot[{XX,YY},{t,0,tmax}]