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


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:

desired result]

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.

  • 3
    $\begingroup$ 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. $\endgroup$
    – Kuba
    Oct 30, 2019 at 9:16
  • $\begingroup$ 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 $\endgroup$
    – isaac5122
    Oct 30, 2019 at 9:21
  • $\begingroup$ Unless I am mistaken, make sure that the initial net momentum of the system is 0, otherwise it will drift. $\endgroup$
    – Kuba
    Oct 30, 2019 at 9:27
  • 1
    $\begingroup$ 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? $\endgroup$
    – yarchik
    Oct 30, 2019 at 16:20
  • 1
    $\begingroup$ "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. $\endgroup$ Nov 29, 2019 at 22:12

1 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.

  • 3
    $\begingroup$ But why do you speak about the three masses? The question clearly states it's a 2-body problem. $\endgroup$
    – yarchik
    Oct 30, 2019 at 16:15
  • $\begingroup$ The two body problem has closed solution. It is not necessary to solve differential equations for this case as the previously shown graphic from the questioning person proved. Only in a more general case for in principle unbound masses the differential equations have be solved. Since the center of mass induces motion on its own in general it can be treated as an additional mass. The is then better solved with the relative intermass distance as a coordinate. The picture already removed misguided in which the two masses orbited their center of mass. $\endgroup$ Nov 1, 2019 at 13:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.