I want to simulate the orbits of planets and moons in our solar-system. They will form a carousel of gravitational bodies that my little spaceships will lift-off from, orbit around, and slingshot past, on their imaginary journeys of adventure and derring-do.

The wikipedia article on orbit modeling is a nice introduction to some of the concepts involved. Newton’s laws form the basis of computing orbits – every body attracts every other body in proportion to the product of their masses and inverse to the square of the distance between them.

If you imagine two bodies (say the Sun and Earth) orbiting in isolation, then an analytic solution (a directly computable equation) to the orbit is possible. It’s a conic section (circle, ellipse or parabola) called a Kepler orbit.

A Kepler orbit, being a nice geometric figure, can be described rather simply with a small set of parameters.

Of course, the Sun and the Earth aren’t the only two bodies involved. Though the Sun’s mass dwarfs all others, planets and moons also pull on each other causing deviations – called perturbations – from the Keplerian orbit, which are evident when you measure carefully enough. There are no analytical solutions to the general problem that considers more than two bodies. Orbit prediction (also called propagation) is done using a variety of methods with the complexity dependent on the level of accuracy needed.

The most direct way that comes to mind – which is what I have used in previous code – is formally called an N-body integrator. Here, if you have N bodies in the simulation you simply compute the force F between every pair of bodies, vector sum the resultant force for each body, find the accelerations, and then use numerical integration to update the velocities and positions. I’ll call this the N-body method. The basic principle is straightforward, but there are important details to do with the numerical integration itself.

NASA’s HORIZONS system (pdf introduction) is the goto place for ephemeris data (ephemeris means “diary” in Latin and Greek and is the historical name for tables of positions of astronomical objects). The whole HORIZONS system has a quaint 90s scientific computer system feel to it. There are three access methods – via web, via email and via telnet. Yes, telnet. I haven’t heard that word in a decade. I was pretty surprised when I found it still existed on my mac.

The web interface is the nod to “modernity”. It allows you to experiment with queries and familiarize yourself with the system. The web interface is hard to scale, and I started to look for a Python based programatic interface to the HORIZONS system. My first hit was callhorizons. Callhorizons seems to have been created with a lot of love but it returns data in a somewhat specialized format. Scanning through the code and poking a little around the batch web interface it turns out it’s quite simple to create a programatic interface to query the system. I’ve put some barebones Python code to query the HORIZONS system here.

So, another method to determine orbits of planets for the simulation would be to get ephemeris data (I had to get used to that term too) from NASA’s Horizons system and store it as a look-up table. In this case the planetary orbits are completely predetermined and given with a high degree of fidelity. I’ll call this the look-up table method.

After I started doing my research I came across the paper: Keplerian Elements for Approximate Positions of the Major Planets. Here the author takes high accuracy data over a certain time period (epoch) from NASA and fits Keplerian orbits to the data – Keplerian orbits are ellipses and are completely defined by a small number of parameters (called elements). One neat thing the author does is to include delta terms that change the Keplerian elements with time – our ellipses precess round slowly – and this more or less captures the perturbation of the planetary orbits due to the presence of other planets. I’ll call this the “keplerian fit” approach.

My inclination is to start with the N-body simulation approach which is what we’ll explore in the next post.