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.
EllipticE[]
is built-in (and quite fast), so you really don't need to pull out an approximation. $\endgroup$