I would like to build a gravity model (with 3D graphics) that would build the orbit of the planets based on orbital elements input and the gravity calculation for each element.

My initial code is live here (Click the 0.005 to turn it into 0.05 and then click on restart time to run the orbits).

The example above has all planets on Solar System (up to Neptune), the Sun and the asteroid Apophis. For clarity, the Sun and planets other than earth have been hidden, but are still considered for gravity calculations.

However, the orbits do not match the correct estimate intersection date. In this article, it says that Earth-Apophis close-encounters will happen on 13 April (meaning that the earth should be in the orbit's intersection point every 13 April, or very close to it).

The data from the program is taken from Nasa JPL.

Start date is 2454441.5, as given in JPL Apophis data. Planets orbital elements have been corrected to this date, with the calculations in the above PDF.

Why the orbits intersection date is wrong? Should I be doing something else? Am I doing something wrong?

The above Apophis link also cites a "Standard Dynamical Model". If anyone knows how to correctly create one, I believe this would answer my question as well.

TLDR: How can I correctly create planetary orbits from kepler elements?

Thank you in advance if anyone can help me!

  • $\begingroup$ Answer to a deleted comment: I'm converting keplerian orbital elements into cartesian coordinates and then solving it numerically (n^2 total gravity calculations, n = number of bodies). What I meant by the date is that Earth should be at the intersect position on April 13 every year. Apophis position only matters on 2029 and 2036 like you said. $\endgroup$
    – hawaii
    Aug 15, 2016 at 15:43
  • $\begingroup$ I don't know how to do that, but for all 8 planets and Apophis, the calculation is enormously complicated. Once you have 3 bodies, simple equations no longer work and you need to run constant corrections. There's no longer a simple formula. (look up 3 body problem). I don't think what you ask can be programmed without a PHD in perturbation theory and pages and pages of code and perhaps a supercomputer. $\endgroup$
    – userLTK
    Aug 15, 2016 at 18:41
  • 1
    $\begingroup$ @userLTK No, straight (numerical) integration of the system should be accurate enough Eg Verlet integration. I'd guess there is some mistake in the initialisation of the system. $\endgroup$
    – James K
    Aug 15, 2016 at 19:32
  • $\begingroup$ Thanks @JamesK I didn't know that. I remember the 3 body problem being unsolvable when I studied math, but I'm guessing there's work-arounds with a stable system. $\endgroup$
    – userLTK
    Aug 16, 2016 at 11:02

1 Answer 1


If you're running your own N-body dynamics, you can get initial Cartesian positions and velocities for each body in your system (including the Sun) by asking JPL HORIZONS for Ephemeris Type: Vectors. I would choose Coordinate Origin: Solar System Barycenter so that the Sun does not drift away from the origin, and Target Body: [Planet] Barycenter (or Sol) unless you are modeling the Earth and Moon as separate bodies.

If you still want to compute initial Cartesian positions and velocities from Keplerian elements, first do a sanity check: Earth should cross the X axis at the equinoxes and the Y axis at the solstices. Second, those approximate elements for the major planets may not be precise enough for your purpose. Instead you might ask HORIZONS for Ephemeris Type: Elements, and Time Span beginning with the epoch of your elements for Apophis. You could also get elements for Apophis in a more convenient epoch this way.

Note that Keplerian elements are subject to dynamic evolution and are only valid for a limited time near their epoch. If you ask HORIZONS for Apophis's orbital elements one week before and after its close approach to Earth, you can see that they are quite different. Of course the latter are also highly uncertain.

For the gory post-Newtonian details of the dynamical model the Apophis prediction article discussed extending, see JPL's notes on DE430 and DE431, part III.

  • $\begingroup$ I tried to follow closely your answer, but could not fix my code to achieve the first sanity check (very nice test indeed, thanks for that). Can you provide a link to an article or book with the calculation/explanation? $\endgroup$
    – hawaii
    Aug 16, 2016 at 3:07
  • $\begingroup$ My sanity check is based on the usual definition of rectangular ecliptic coordinates. $\endgroup$
    – Mike G
    Aug 16, 2016 at 21:11
  • $\begingroup$ Your answer helped me a lot to debug my code. It is working now. Strangely, the orbital elements (planets) on the NASA page seems to be acting weird. I'm implementing a way to query horizons through telnet and respond to the browser as a restful application. After that I believe I will be able to answer my own question on what was going wrong, since I will have all the tools I need to check everything. Thanks again! $\endgroup$
    – hawaii
    Aug 16, 2016 at 22:14
  • $\begingroup$ The approximate elements for the planets are approximately correct. Earth may have wildly different $\Omega$ and $\omega$ due to near-zero $I$, but the longitude of perihelion should be about the same. Good luck! $\endgroup$
    – Mike G
    Aug 17, 2016 at 13:45

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .