The Classic 3-Body Problem in a Cataclysmic Variable Star System

friction=0.06;
Vx0=-.02;
Vy0=0;
q=.4;
inclination=80;

μ=((M1 q)/(M1+q M1));a=1;M1=1;G=1;Porb=(Sqrt[((4 (π^2) (a^3))/(G (1+q)))]);Ωorb=((2 π)/Porb);M2=q M1;Φ[x_,y_,z_]:=-((G M1)/(Sqrt[(((x+(q/(1+q))))^2)+(y^2)+(z^2)]))-((G M1 q)/(Sqrt[(((x-a+(q/(1+q))))^2)+(y^2)+(z^2)]))-(1/2) ((Ωorb)^2) ((((x-μ+((q a)/(1+q))))^2)+(y^2));deriv=(∂_x)Φ[x,0,0];posAtL1=FindRoot[deriv==0,{x,0.3`}];ΦAtL1=Φ[x/.posAtL1,0,0];posAtL2=FindRoot[deriv==0,{x,-1}];ΦAtL2=Φ[x/.posAtL2,0,0];posAtL3=FindRoot[deriv==0,{x,1.5`}];ΦAtL3=Φ[x/.posAtL3,0,0];posAtL4=FindRoot[(∂_y)Φ[x/.posAtL1,y,0]==0,{y,1}];ΦAtL4=Φ[x/.posAtL1,y/.posAtL4,0];μ=mu/.Solve[q==(mu/(1-mu)),mu][[1]];grL1=ContourPlot3D[Φ[x,y,z],{x,x/.posAtL1,x/.posAtL3},{y,-0.6`,0.6`},{z,-0.6`,0.6`},Contours->{ΦAtL1},ContourStyle->{Red,Specularity[White,3]},BoxRatios->Automatic,Mesh->None];solution=NDSolve[{x''[t]==2 y'[t]+x[t]-(((1-μ) (x[t]+μ))/((((((x[t]+μ))^2)+((y[t])^2)))^(3/2)))-((μ (x[t]-1+μ))/((((((x[t]-1+μ))^2)+((y[t])^2)))^(3/2)))-friction x'[t],y''[t]==-2 x'[t]+y[t]-(((1-μ) y[t])/((((((x[t]+μ))^2)+((y[t])^2)))^(3/2)))-((μ y[t])/((((((x[t]-1+μ))^2)+((y[t])^2)))^(3/2)))-friction y'[t],y[0]==0,x[0]==Evaluate[x/.posAtL1],y'[0]==Vy0,x'[0]==Vx0},{x[t],y[t]},{t,0,13}];

lst=Table[Show[ParametricPlot3D[Evaluate[{x[t],y[t],0}/.solution],{t,0,orbtime},PlotStyle->{Red}],grL1,Graphics3D[{PointSize[0.02`],White,Point[{-μ,0,0}],Text["M1",{-μ,0,-0.1`}],Text["M2",{1-μ,0,-0.1`}],Text["CM",{0,0,-0.1`}],Text["+",{0,0,0}]},Axes->False],AxesLabel->{"x","y","z"},PlotRange->{{x/.posAtL2,x/.posAtL3},{-0.8`,0.8`},{-0.6`,0.6`}},Lighting->{{"Ambient",GrayLevel[.3]},{"Point",White,{-μ,0,0}}},Background->Black,SphericalRegion->True,Axes->False,Boxed->False,PlotLabel->Style[Column[{Row[{"q = ",ToString[q],"   friction = ",ToString[friction]}],Row[{Subscript["V","x"],"0 = ",ToString[Vx0],"   ",Subscript["V","y"],"0 = ",ToString[Vy0]}],Row[{"orbits = ",ToString[orbtime]}]}],14,White],ImageSize->{400,300},ViewPoint->{2 Cos[-2 π orbtime] Sin[((inclination π)/(180.`))],2 Sin[-2 π orbtime] Sin[((inclination π)/(180.`))],Cos[((inclination π)/(180.`))]},ViewAngle->π/12],{orbtime,0.1`,3,0.01`}];

SetDirectory[NotebookDirectory[]];

Export["threebody3danim.mp4",lst]