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

<< Graphics`

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

Table[
μ = 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/(1 + q) a))^2) + 
                      (y^2));
deriv = D[Φ[x, 0, 0], x];
posAtL1 = 
              FindRoot[deriv == 0, {x,  .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[
                D[Φ[x /. posAtL1, y, 0], y] == 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, \(- .6\),  .6}, {z, \(- .6\),  .6}, 
				PlotPoints -> {7, 7}, Contours -> {ΦAtL1}, 
                ContourStyle -> {EdgeForm[], 
                    SurfaceColor[Red, White, 3]}, 
                BoxRatios -> Automatic, 
                DisplayFunction -> Identity]\);
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}];
\(Show[
              ParametricPlot3D[
                Evaluate[{x[t], y[t], 0, {RGBColor[1, 0, 0]}} /. 
                    solution], {t, 0, orbtime}, 
                DisplayFunction -> Identity], grL1, 
              Graphics3D[{PointSize[0.02], White, Point[{\(-μ\), 0, 0}], 
                  Text["M1", {\(-μ\), 0, \(- .1\)}], 
                  Text["M2", {1 - μ, 0, \(- .1\)}], 
                  Text["CM", {0, 0, \(- .1\)}], 
                  Text["+", {0, 0, 0}]}, Axes -> False, 
                DisplayFunction -> Identity], 
              TextStyle -> {FontFamily -> "Arial", 
                  FontWeight -> "Bold"}, 
              AxesLabel -> {"x", "y", "z"}, 
              PlotRange -> {{x /. posAtL2, 
                    x /. posAtL3}, {\(- .8\),  .8}, {\(- .6\),  .6}}, 
              LightSources -> {{{Cos[\(-2\)\ π\ orbtime - π/
                              2]\ Sin[inclination\ π/180. ], \(-Cos[
                          inclination\ π/180. ]\), 
                      Sin[\(-2\)\ π\ orbtime - π/2]\ Sin[
                          inclination\ π/180. ]}, White}}, 
              AmbientLight -> RGBColor[ .3,  .3,  .3], 
              Background -> Black, SphericalRegion -> True, 
              Axes -> False, Boxed -> False, 
              PlotLabel -> 
                StyleForm["q=" <> ToString[q] <> "   friction=" <> 
                    ToString[friction] <> \*""\n\!\(Vx\_0\)="" <> 
                    ToString[Vx0] <> \*""   \!\(Vy\_0\)="" <> 
                    ToString[Vy0] <> "
orbits=" <> ToString[orbtime], 
                  FontSize -> 14], ImageSize -> 600, 
              ViewPoint -> {2\ Cos[\(-2\)\ π\ orbtime]\ Sin[
                      inclination\ π/180. ], 
                  2\ Sin[\(-2\)\ π\ orbtime]\ Sin[
                      inclination\ π/180. ], 
                  Cos[inclination\ π/180. ]}, 
              DisplayFunction -> $DisplayFunction]\),{orbtime, 0.1, 3,  .1}];