Gravitational Potential of a Cataclysmic Variable Star System

First we define some initialization constants

q =  .4;
μ = 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;

Next, we define the potential

Φ[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))

Here, we determine the location of several critical points, namely the Lagrange points

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];

Gas Trajectories and the 3 Body Problem

Here we define the initial velocities of the incoming gas particles and a coefficient of linear friction

Vx0 = -0.01;
Vy0 = 0.01;
friction = 0.08;

Next, we construct the shape of the secondary star, determine the path of the incoming gas stream, and combine the results

μ = 
        mu /. Solve[q == mu/(1 - mu), 
              mu][[1]];
\(grL1 = 
        First[Graphics[
            ContourPlot[Φ[x, y, 0], {x, x /. posAtL1, 
                1.5}, {y, \(- .6\),  .6}, 
				PlotPoints -> 50, 
              Contours -> {ΦAtL1}, 
              ContourShading -> False, AspectRatio -> 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[
      Graphics[{\(ParametricPlot[{x[t], y[t]} /. solution, {t, 0, 10}, 
              PlotStyle -> {RGBColor[1, 0, 0]}, 
              DisplayFunction -> Identity]\)[\([1]\)], PointSize[0.02], 
          Point[{\(-μ\), 0}], Text["M1", {\(-μ\), \(- .1\)}], 
          Text["M2", {1 - μ, \(- .1\)}], 
          Text[StyleForm[\*""\!\(Vx\_0\)="" <> ToString[Vx0], 
              FontSize -> 14], { .5,  .55}], 
          Text[StyleForm[\*""\!\(Vy\_0\)=+"" <> ToString[Vy0], 
              FontSize -> 14], { .5,  .45}], 
          Text[StyleForm["friction=" <> ToString[friction], 
              FontSize -> 14], {\(- .5\),  .55}], 
          Text["CM", {0, \(- .1\)}], Text["+", {0, 0}], grL1}, 
        Frame -> True, 
        PlotRange -> {{\(-1.2\), 1.2}, {\(- .6\),  .6}}, 
        AspectRatio -> 1, ImageSize -> 400, 
        PlotLabel -> 
          StyleForm["q=" <> ToString[q], FontSize -> 14], 
        Axes -> False], 
      TextStyle -> {FontFamily -> "Times", 
          FontWeight -> "Bold"}, FrameLabel -> {"x", "y"}, 
      ImageSize -> 800, AspectRatio -> Automatic]\)

\(\(secondary = 
        First[Graphics[
              ContourPlot[Φ[x, y, 0], {x, x /. posAtL1, 
                  1.5}, {y, \(- .6\),  .6}, 
				PlotPoints -> 50, 
                Contours -> {ΦAtL1}, 
                ContourStyle -> {Red}, ContourShading -> False, 
                AspectRatio -> Automatic, 
                DisplayFunction -> Identity]]] /. 
          Line[x__] -> Polygon[x];\)\)

\(\(anim = 
        Table[Show[
            Graphics[{White, Text["M1", {\(-μ\), \(- .1\)}], 
                Text["M2", {1 - μ, \(- .1\)}], 
                Text[StyleForm[\*""\!\(Vx\_0\)="" <> ToString[Vx0], 
                    FontSize -> 14], { .5,  .45}], 
                Text[StyleForm[\*""\!\(Vy\_0\)=+"" <> ToString[Vy0], 
                    FontSize -> 14], { .5,  .35}], 
                Text[StyleForm["friction=" <> ToString[friction], 
                    FontSize -> 14], {\(- .5\),  .45}], 
                Text["CM", {0, \(- .1\)}], Text["+", {0, 0}], 
                Point[{\(-q\)/\((1 + q)\), 0}], secondary, Red, 
                Table[Point[{x[t], y[t]} /. solution[\([1]\)]], {t, 0, 
                    time,  .01}]}], Frame -> True, 
            PlotRange -> {{\(-1.2\), 1.2}, {\(- .6\),  .6}}, 
            AspectRatio -> Automatic, ImageSize -> 400, 
            Background -> Black, 
            TextStyle -> {FontFamily -> "Arial", 
                FontWeight -> "Bold", 
                FontColor -> White}], {time, 0, 13,  .1}];\)\)

Here we set the working directory to the notebooks directory and export the animation as an animated GIF.

SetDirectory[
        ToFileName[("FileName" /. 
                NotebookInformation[
                  EvaluationNotebook[]])[[1]]]];

Export["3bodyanim.gif", anim, "GIF", 
      ConversionOptions -> {"Loop" -> True}]

Isosurfaces

\(Plot[Φ[x, 0, 0], {x, \(-1.5\), 2}, 
	
      PlotRange -> {0, \(-15\)}, 
	
      FrameLabel -> {"x/a", "Potential"}, Frame -> True, 
      TextStyle -> {FontWeight -> "Bold", 
          FontSlant -> "Italic", FontFamily -> "Times"}, 
      ImageSize -> 550]\)

\(Plot3D[Φ[x, y, 0], {x, \(-1.5\), 2}, {y, \(-1.5\), 
        1.5}, 
	PlotRange -> {0, \(-5\)}, 
	Boxed -> False, \ 
	
      Axes -> False, 
	ViewPoint -> {2.841, \ 0.144, \ 1.833}, 
	
      PlotPoints -> 50, ImageSize -> 550, 
      AmbientLight -> RGBColor[0.43, 0.35, 0.5]]\)

<< Graphics`

\(grL1 = 
      Show[Insert[
          ContourPlot3D[Φ[x, y, z], {x, x /. posAtL3, 
              x /. posAtL2}, {y, \(- .6\),  .6}, {z, \(- .6\),  .6}, 
				PlotPoints -> {11, 11}, 
				Boxed -> False, 
            Contours -> {ΦAtL1}, 
            ContourStyle -> {{SurfaceColor[Orange, White, 
                    20]}}, 
				ViewPoint -> {\(-0.011\), \(-3.199\), 1.104}, Axes -> False, AmbientLight -> Black, 
            DisplayFunction -> Identity], EdgeForm[], {1, 1}], 
        DisplayFunction -> $DisplayFunction]\)

\(grL2 = 
      Show[Insert[
          ContourPlot3D[Φ[x, y, z], {x, \(-1.6\), 1.6}, {y, 0, 
              1.6}, {z, \(- .7\),  .7}, 
				PlotPoints -> {11, 
                11}, 
				Boxed -> 
              False, 
				Contours -> {ΦAtL2}, 
            ContourStyle -> {{SurfaceColor[Blue, White, 
                    20]}}, 
				ViewPoint -> {\(-0.011\), \(-3.199\), 1.104}, Axes -> False, AmbientLight -> Black, 
            DisplayFunction -> Identity], EdgeForm[], {1, 1}], 
        DisplayFunction -> $DisplayFunction]\)

\(grL3 = 
      Show[Insert[
          ContourPlot3D[Φ[x, y, z], {x, \(-1.6\), 1.6}, {y, 0, 
              1.6}, {z, \(- .7\),  .7}, 
				PlotPoints -> {5, 
                5}, 
				Boxed -> 
              False, 
				Contours -> {ΦAtL3}, 
            ContourStyle -> {{SurfaceColor[RGBColor[0.08, 0.85, 1], 
                    White, 20]}}, 
				ViewPoint -> {\(-0.011\), \(-3.199\), 1.104}, Axes -> False, AmbientLight -> Black, 
            DisplayFunction -> Identity], EdgeForm[], {1, 1}], 
        DisplayFunction -> $DisplayFunction]\)

\(grL4 = 
      Show[Insert[
          ContourPlot3D[Φ[x, y, z], {x, \(-1.6\), 1.6}, {y, 0, 
              1.6}, {z, \(- .7\),  .7}, 
				PlotPoints -> {5, 
                5}, 
				Boxed -> 
              False, 
				Contours -> {ΦAtL4 -  .03}, 
            ContourStyle -> {{SurfaceColor[Purple, White, 
                    20]}}, 
				ViewPoint -> {\(-0.011\), \(-3.199\), 1.104}, Axes -> False, AmbientLight -> Black, 
            DisplayFunction -> Identity], EdgeForm[], {1, 1}], 
        DisplayFunction -> $DisplayFunction]\)

myout = 
      Show[grL1, grL2, grL3, grL4, ImageSize -> 800, 
        LightSources -> {{{1, 1, 1}, White}, {{-1, 1, 1}, 
              White}}]