2 body orbits around a barycenter

Question

I'm attempting to plot a barycentric orbit of the sun and Jupiter.

Initial Conditions

T = QuantityMagnitude[UnitConvert[PlanetData["Jupiter","OrbitPeriod"]], 
"days"](*period of jupiter*)
ψ = {180, 0}(*True Anomoly*)
Ecc = {0.04839266};(*Eccentricity of jupiter*)
a = {5.2033630096869589997`8.99956592252068};(*Semi major and minor axis of jupiter*)
b = {5.1972666917898543459`8.989499285086092};
m = {1.`,0.0009546133303706552`};(*mass of sun and jupiter in solar masses*)
G = 0.00029589743849552926`;(*gravitational constant*)
μ = G*m;(*standard gravitational parameter of the sun and jupiter*)
Subscript[x, cm] = (a*m[[2]])/(m[[1]])(*centre of mass of the system in terms of semi major axis, this is also the semi major axis of the sun from barycenter*)
Subscript[y, cm] = (b*m[[2]])/(m[[1]])(*centre of mass of the system in terms of semi minor axis, this is also the semi minor axis of the sun from barycenter*)
Subscript[Ecc, sun] = Sqrt[1 - (b/a)^2](*eccentricity of the sun*)
Subscript[a, j] =a - Subscript[x,cm](*updated semi major axis of jupiter: this shows the semi major axis from the barycentre instead of from the centre of mass of the sun*)
Subscript[cm, a] = Flatten[{Subscript[x, cm], Subscript[a,j]}](* semi major axis of sun and jupiter from barycenter*)
Subscript[cm, Ecc] = Flatten[{Subscript[Ecc, sun],Ecc}](* Eccentricity of sun and jupiter*)

Calculating the orbital position and orbital velocity of the sun and Jupiter around barycenter

The orbital position in terms of the true anomaly is given by

$$r=\frac{a(1-e^2)}{1+e\cos[\psi]}$$

Where $a$ is semi major axis, e is the eccentricity and $\psi$ is the true anomaly

The respective orbital velocity is given as $ $v=\sqrt{\mu(\frac{2}{r}-\frac{1}{a})}$$

When calculating the velocity of the sun around barycenter, I assumed that it had the same period of Jupiter. Therefore I approximated the elliptical path taken using the semimajor and minor axis. A handy calculator can be seen here.

If the period of Jupiter is T=4332.8201 and the approximate path of the sun is d= 0.03119155667 Astronomical units(AU) Then the orbital velocity of the sun in terms of the period of Jupiter is: v=d/t v=7.200267006001846*10^-6 AU/day

r = Table[Subscript[cm, a][[i]] (1 - Subscript[cm, Ecc][[i]]^2)/(1 + 
  Subscript[cm, Ecc][[i]] Cos[ψ[[i]] Degree]), {i,2}] (*orbital position, From this, the sun should be at the left side of the barycentre and jupiter should be at the right hand side*)
rx = Table[r[[i]] Cos[ψ[[i]] Degree], {i, 2}](*x component of position*);
ry = Table[r[[i]] Sin[ψ[[i]] Degree], {i, 2}](*y component of position*);
Subscript[v, jupiter] = Sqrt[μ[[1]] (2/r[[2]] -1/a[[1]])](*orbital velocity of jupiter at respective true anomoly*)
v = {-7.200267006001846`*^-6,0.007922399185159456`}(*updated velocity such that the sun orbits the path with respect to the period of jupiter*)
vx = -v Sin[ψ Degree](*x component of velocity*);
vy = v Cos[ψ Degree](*y component of velocity*);

Solving equations of motion and plotting

eq = {Table[x[i]''[t] == Sum[If[j == i,0, (-μ[[j]] (x[i][t] - 
        x[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - 
          y[j][t])^2)^(3/2)], {j, 2}], {i, 2}],Table[y[i]''[t] == 
 Sum[If[j == i,0, (-μ[[j]] (y[i][t] - y[j][t]))/((x[i][t] - x[j][t])^2 + (y[i][t] - y[j][t])^2)^(3/2)], {j, 2}], {i, 2}]};
var = Join[Table[x[i], {i, 2}], Table[y[i], {i, 2}]];
orb = NDSolve[{eq, Table[x[i][0] == r[[i]], {i, 2}], 
Table[y[i][0] == 0, {i, 2}], Table[x[i]'[0] == 0, {i, 2}], 
Table[y[i]'[0] == v[[i]], {i, 2}]}, var, {t, 90600}];
plot2D = Show[Table[ParametricPlot[Evaluate[{x[i][t], y[i][t]} /. orb], {t, 0, T}, 
 PlotStyle -> Blue, PlotRange -> 6], {i, 2}]];
Animate[Show[plot2D,Graphics[Table[{Red, PointSize[0.02], 
 Point[{x[i][t], y[i][t]} /. orb]}, {i, 2}]]], {t, T}, AnimationRate -> 50, AnimationRunning -> False]

The problem

When plotting this at range of 6(AU) I receive this.

and all may appear correct. However upon closer inspection by changing PlotRange to 0.01 and changing values of t to a range e.g t,0,10000

I receive this.

I'm a little confused to why this happening, I'm trying to achieve something like this:

where the sun orbits the barycenter in an epicyclodic path . No matter what I try I cant stop the sun from drifting up!! I hope I've provided enough information! Can anyone help solve this?

It should be noted that this question is an extension of my previous: gravitational two body problem for the orbit of the sun and jupiter around their barycenter

In my opinion I did not provide enough information hence the new post. Apologies if this is a duplicate thread.

Answer 1

Gravitation is often regarded as conservative or to pose several invariants. This system of three gravitational masses do have energy conservation. So You start from the very right thoughts. But Your physical power is not large enough to digg deep into the system. So Kuba is right. The three masses should in general move. They move each.

The solution to control the rest of one of them is to add condition to the system of differential equation. It is on the other hand not trivial to find the conserved quantities in the system. The first one already implemented, the three masses move in one common plane for all time, might make it even harder. The condition all are treated as mass point with no internal degrees of freedom is unnatural as well. For a first approximation try the most general physical quantities in the system. For example the center of mass, the moment of rotation relativ to the center of mass. Control the conservation of energy of the degree of energy dissipation. The conservation of energy might be given for the whole system in motion, the moving system each or pairs of them. Since the main idea is classical gravitational interaction in simplified differential equation their might be plenty of space to variate and calculate.

In the publications from Wolfram for Mathematica 12 there is an example for three body motion: NBodySimulation. This shows up the newest trend dealing with such kind of physical problems in Mathematica. So have fun.