# 2 body orbits around a barycenter

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.20336300968695899978.99956592252068};(*Semi major and minor axis of jupiter*)
b = {5.19726669178985434598.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[])/(m[])(*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[])/(m[])(*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[μ[] (2/r[] -1/a[])](*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] == r[[i]], {i, 2}],
Table[y[i] == 0, {i, 2}], Table[x[i]' == 0, {i, 2}],
Table[y[i]' == 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'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.

• Is this related to your recent question: 208467? You've received feedback there and left it without reply so that does not encourage users to help next time.
– Kuba
Oct 30, 2019 at 9:16
• Mea culpa, I wasnt aware my feedback had not posted, thank you for highlighting this. This is in reference to the previous question however i felt as though there was not enough detail and explanation to how i derived the constants and arrived at the problem. I hope this clarifies the situation Oct 30, 2019 at 9:21
• Unless I am mistaken, make sure that the initial net momentum of the system is 0, otherwise it will drift.
– Kuba
Oct 30, 2019 at 9:27
• I have a problem with the way you organize your calculation. First you have an initialization part. Very good....From now on you should only use these parameters. However, in the second block one suddenly sees "updated velocity such that the sun orbits the path with respect to the period of jupiter". What is this and why it is injected there? Why don't you stick to your parameters defined above? Oct 30, 2019 at 16:20
• "I approximated the elliptical path taken using the semimajor and minor axis." - EllipticE[]` is built-in (and quite fast), so you really don't need to pull out an approximation. Nov 29, 2019 at 22:12