I have developed an enthusiasm (aka weird obsession) for celestial mechanics and developed several games and the Unity Asset Gravity Engine . In developing Gravity Engine I learned a great deal about high pedigree N-body simulations, but in some cases (e.g. a model of the solar system) it is not really necessary to simulate the system but rather just evolve it in the correct way. In this case Gravity Engine offers the option to use Kepler’s equation and move bodies in their elliptical orbit as a function of time. This uses far less CPU than doing the 9*8 mutual gravitational interactions (10*9, if you add in the “dwarf planet” Pluto).
[If you don’t see an animation, click the post title to see ONLY this post. There is WordPress JS bug when this animation and a YouTube link are present]
Creating code to move a planet in an elliptical orbit with the correct velocity is surprisingly tricky. This is one of those cases where you might expect you could just grab a formula from Wikipedia and bang out some code. This bias comes from all the examples in physics class where the goal is to “find/derive the formula” and get a tidy equation.
If you dig around on physics pages for an equation of an elliptical orbit you will generally encounter the equation for the shape of the orbit with eccentricity e, and semi-major axis, a:
I have attached the subscript F to the angle to indicate this is the angle from the focus of the ellipse between the position of the body and the long axis of the ellipse. Historically, this angle is named the True Anomaly.
Where is time in this equation?
Nowhere. This equation doesn’t tell us anything about time.
To get an equation for how the object moves as a function of time, we’ll need Kepler’s equation. Kepler constructed his equation without Calculus (which came along about 60 years after Kepler did this work) using geometric arguments and the assumption that an objects speed in an elliptical orbit was inversely proportional to it’s distance from the focus. Kepler’s equation is:
Here M is the position of a body moving in a uniform circle at a constant rate (that we will relate to time) and E is the angle to the position on the ellipse from the origin, called the eccentric anomaly. Here a picture will help:
The eccentric anomaly 
If we have a specific time we want a position for we need to convert this into a value of M. This is done by dividing the time by the time per orbit, T. Kepler can also help us here with his third law that relates the size and eccentricity of the orbit to the mass of the bodies:
where m is the combined mass of the central object and orbiting body. (Kepler did not know the proportionality constant was the mass, that came later).
Given M, Solve for E
Ok, we’ll just isolate E…hmmm. E appears by itself and inside the sine function. That sinks our chance of getting a tidy mathematical formula. This equation is legendary in Mathematics, since it is an early example of a transcendental equation with an import application. It has been studied extensively and the approaches are well summarized in the book “Solving Kepler’s Equation Over Three Centuries” by Peter Colwell.
There are some series approximations, but they are not valid for all eccentricities. The most common approach is to iterate the equation until we converge on a value that is “good enough”.
The “recipe” for tying this all together is:
- Determine the orbital period T
- For the time t we’re interested in, divide by the orbital period and use the remainder to find the angle if the body were moving in a uniform circle, M.
- Using iteration solve Kepler’s equation and find E
- Using E, determine
- Find the corresponding r using the orbit equation.
In Javascript this becomes:
var orbitPeriod = 2.0 * Math.PI * Math.sqrt(a*a*a/(m*m)); // G=1
function orbitBody() {
// hide last position drawn
context.clearRect( last_x -10, last_y - 10, 20, 20);
// 1) find the relative time in the orbit and convert to Radians
var M = 2.0 * Math.PI * time/orbitPeriod;
// 2) Seed with mean anomaly and solve Kepler's eqn for E
var u = M; // seed with mean anomoly
var u_next = 0;
var loopCount = 0;
// iterate until within 10-6
while(loopCount++ < LOOP_LIMIT) {
// this should always converge in a small number of iterations - but be paranoid
u_next = u + (M - (u - e * Math.sin(u)))/(1 - e * Math.cos(u));
if (Math.abs(u_next - u) < 1E-6)
break;
u = u_next;
}
// 2) eccentric anomoly is angle from center of ellipse, not focus (where centerObject is). Convert
// to true anomoly, f - the angle measured from the focus. (see Fig 3.2 in Gravity)
var cos_f = (Math.cos(u) - e)/(1 - e * Math.cos(u));
var sin_f = (Math.sqrt(1 - e*e) * Math.sin (u))/(1 - e * Math.cos(u));
var r = a * (1 - e*e)/(1 + e * cos_f);
time = time + 1;
// animate
last_x = focus_x + r*cos_f;
last_y = focus_y + r*sin_f;
drawBody( r* cos_f, r*sin_f, "blue");
setTimeout(orbitBody, 30);
}
I have left out some of the init code for clarity. If you view the source for this page you can find all this.
(Eccentric anomaly image created by CheCheDaWaff/Creative Commons License).
are needed to describe a point. If we think of the typical donut representation sitting on the x, y plane then 
is the angle around the z-axis (say from the line x=0) and 
is the angle around the cross section of the torus at a given 
and the radius of the cross section is 
.
is the same at the origin as it is at any other point. On curved surfaces things are different. To make this clear lets consider a simpler 2D surface, the sphere, with longitude 
and latitude 
on a curved surface come from a branch of mathematics known as differential geometry. “Differential” since it studies the behaviour of surfaces that are smooth enough that you can do calculus, i.e. no sudden steps or cusps. One of the new ideas that enters with differential geometry is that while you can do calculus at each point, you cannot directly compare the results of these calculations at different points. This is because each point has its own tangent plane (more generally tangent space). As an object slides along the 2D surface the local co-ordinates of the object change relative to the surface co-ordinates. Returning to our sphere example, if you initially were heading west from New York when your great circle crosses the equator you are no longer facing west BUT you did not turn, you were just following your geodesic path. Exactly what the geodesic path is and how it turns as you go is something that can be calculated for a surface.
. This was a good introduction to OpenGL ES and I went between it and the 
. I was able to get some reasonable 3D stuff going. I did find the C++ cross-platform framework in the iPhone book a bit unsatisfying. After all I was also trying to get some chops in Objective-C. However, I had running code and started to work on the touch interface and game flow.
seriously. This had repercussions in N-body. Since I could not have clean code “at the office” I tried extra hard to keep N-body clean. This resulted in some significant refactoring in N-body.
a big help in getting into the flow of Objective-C coding. Another big bonus was a colleague who had a long history in MacOS development and was generous with his time. I also used this opportunity to go all the way to the bottom in OpenGL. Since I had no strict self-imposed deadline I messed around with how vertices were generated, how the lighting in the shaders work, normal maps for texture simulation, how to map strings into textures into transparent objects, and on and on. (Once iOS6 came along the Xcode OpenGL debugging became much easier. That feature in Xcode is hugely powerful).
and became entranced with the weird tadpole and horseshoe orbits that occur as bodies oscillate about Lagrange points. I really wanted to show these orbits in the co-rotating frame and find a way to interactively generate some intuition about how variations in starting conditions affect the orbit. This added more time to the development. Then I stumbled into the choreography orbits. These were amazing but didn’t quite hang in there with the integrator I had in my physics engine. This sent me on a detour to learn more about choreographies and change the physics engine to do fourth-order Runge-Kutta evolution.