Programming a threebody problem in JavaScript
Please use a newer browser to see the simulation.
In this tutorial we will program motion of three bodies in HTML and JavaScript. This work is built upon the twobody simulation code. I want to say huge thanks to Dr Rosemary Mardling, who taught me astrophysics in Monash University. This work is based on Rosmary’s ideas and code. As always, feel free to check the full source code, and use it for any purpose.
Threebody simulations
The little buttons under the slider run the following simulations.
Figure eight
This is a stable threebody system discovered by Cris Moore [5]. The system remains stable even if we change the masses of all bodies a little bit, to 0.99 for example. Just for fun, try increasing the speed of this animation by clicking the clock icon and moving the slider. At certain speeds you will see weird stroboscopic effects. Be careful, this can make you dizzy.
Sun, Earth and Jupiter
This simulation uses true masses, velocities and distances of the Sun, Earth and Jupiter. We can measure one period of the Earth’s orbit in the simulation to be around one second (may depend on computer speed and refresh rate of the monitor). This means that the simulation is working correctly, because it is run at one year per second.
Lagrange point L5
Here the Earth is located near a special point in space called the SunJupiter L5 Lagrange point. Notice that the radius of the Earth’s orbit is smaller than that of Jupiter initially. Planets that are closer to the Sun have shorter orbital periods. Therefore, normally, the Earth would overtake Jupiter. We can check this by decreasing Jupiter’s mass and clicking the Reload button on the bottom right of the simulation screen. However, the combined gravity from Jupiter and the Sun traps the Earth, and it is destined to remain at L5 point behind Jupiter.
Kepler16
This is a simulation of a binary star system that also has a planet with a mass of 1/3 of Jupiter’s. Both binary stars are smaller than the Sun. The orbital plane of the planet is located edgeon to us. Consequently, when the planet moves in front of the stars, it blocks some of their light. This is precisely how this planet was discovered [6]. The simulation tells us that the system appears to be stable, at least in the short term.
Chaotic
This is an example of how an orderly system can quickly become erratic. When the simulation is running on the computer, small calculation and rounding errors accumulate over time and make the system unpredictable. The simulation is so sensitive to small changes that it may even look different when you run it with the same settings but in different browsers.
The coordinate system
Let’s begin implementing our simulation. We will be using a coordinate system with its origin located at the center of mass of the three bodies, as shown in Fig. 1. We use this origin because, as it turns out, the center of mass of any number of bodies remains motionless or moves with constant velocity.
Figure 1: A coordinate system with the origin at the common center of mass of the three bodies.
The equations of motion
The main equation for this simulation is the Newton’s equation of universal gravitation:
Combining this equation with Newton’s second law F = m a
, we can derive three equations of motion:
Here the double dots above the variables mean second time derivatives. The vector
points from the Sun to the Earth. Equation 2 also includes the magnitudes of the vectors, which can be calculated as follows, for the case of the vector r12
:
Here x1, y1
and x2, y2
are the coordinates of the Earth and the Sun respectively.
Writing equations of motions in x
and y
Equation 2 contains three equations of motion. In order to use the equations in our program we need to write each equation in terms of x
and y
coordinates. This gives six equations
Note the we got rid of the masses on the lefthand sides by dividing both sides by them. The equations look intimidating. However, these are not meant to be used by humans. It is much easier to write these equations in code, since they contain a lot of repetitions. But before we do this, we need to convert the equation into a computer readable form.
Converting the equations to first order
Equation 5 is a system of six second order differential equations. In order to be solved numerically, we need to convert the equations into twelve first order ones. This reduction of order is done by first introducing twelve new variables u
:
Now we can use the new variables and turn Equation 5 into the following system of firstorder differential equations that describe the motion of the three bodies:
Again, this looks horrific, but in the program this turns into a nice loop that takes care of repetitions.
Programming the equations of motion
We are finally done with algebra and will turn our equations into code. The first thing we need to do is to declare an array that will store that values of those twelve u
variables:
// Current state of the system
var state = {
// State variables used in the differential equations
// First two elements are x and y positions,
// and second two are x and y components of velocity
// repeated for three bodies.
u: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
}
Next, as promised, we turn the long Equation 7 into a much shorter code:
// Calculate the derivatives of the system of ODEs that describe equation of motion of the bodies
function derivative() {
var du = new Array(initialConditions.bodies * 4);
// Loop through the bodies
for (var iBody = 0; iBody < initialConditions.bodies; iBody++) {
// Starting index for current body in the u array
var bodyStart = iBody * 4;
du[bodyStart + 0] = state.u[bodyStart + 0 + 2]; // Velocity x
du[bodyStart + 1] = state.u[bodyStart + 0 + 3]; // Velocity y
du[bodyStart + 2] = acceleration(iBody, 0); // Acceleration x
du[bodyStart + 3] = acceleration(iBody, 1); // Acceleration y
}
return du;
}
The function derivative
calculates the derivatives of the twelve u
variables using Equation 7. For simplicity, the righthand sides of equation that calculate accelerations, are coded as a separate function:
// Calculates the acceleration of the body 'iFromBody'
// due to gravity from other bodies,
// using Newton's law of gravitation.
// iFromBody: the index of body. 0 is first body, 1 is second body.
// coordinate: 0 for x coordinate, 1 for y coordinate
function acceleration(iFromBody, coordinate) {
var result = 0;
// Starting index for the body in the u array
var iFromBodyStart = iFromBody * 4;
// Loop through the bodies
for (var iToBody = 0; iToBody < initialConditions.bodies; iToBody++) {
if (iFromBody === iToBody) { continue; }
// Starting index for the body in the u array
var iToBodyStart = iToBody * 4;
// Distance between the two bodies
var distanceX = state.u[iToBodyStart + 0] 
state.u[iFromBodyStart + 0];
var distanceY = state.u[iToBodyStart + 1] 
state.u[iFromBodyStart + 1];
var distance = Math.sqrt(Math.pow(distanceX, 2) +
Math.pow(distanceY, 2));
result += constants.gravitationalConstant *
initialConditions.masses[iToBody] *
(state.u[iToBodyStart + coordinate] 
state.u[iFromBodyStart + coordinate]) /
(Math.pow(distance, 3));
}
return result;
}
There are a lot of arrays and indices in this function, and I’ve made a lot of mistakes when I first coded this. But eventually, all bugs were found, removed and carefully released into the wild.
Solving the equations of motion numerically
Next, we need to solve the Equation 7 numerically. In order to do this, we use a popular numerical method called RungeKutta. At each frame of our animation, we call the rungeKutta.calculate
function and it updates our variables u
with new positions and velocities of the three bodies.
// The main function that is called on every animation frame.
// It calculates and updates the current positions of the bodies
function updatePosition(timestep) {
rungeKutta.calculate(timestep, state.u, derivative);
}
Before the first run, however, we need to choose some initial positions and velocities of the Sun, Earth and Jupiter and store them in the u
array.
Drawing the bodies on screen
At each frame of the animation, we calculate new positions of the three bodies in the u
array. These position are turned into coordinates of the bodies on computer screen. Here we use the same drawing techniques that we discussed in the harmonic oscillator tutorial.
📚📖✨
If you are into science fiction, The ThreeBody Problem novel by Liu Cixin and the other two books of the trilogy could be worth your time.
Credits

This work is based on code and lectures by Dr Rosemary Mardling from Monash University.

“The Blue Marble” image: NASA/Apollo 17 crew; taken by either Harrison Schmitt or Ron Evans, source, source.

“The Sun photographed at 304 angstroms” image: NASA/SDO (AIA), source, source.

“Jupiter’s South Pole” image: NASA/JPLCaltech/SwRI/MSSS/Betsy Asher Hall/Gervasio Robles, source.

Figure eight orbit: Moore, C. 1993, Phys. Rev. Lett., 70, 3675.

Kepler16 system: Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602.
References
 The complete source code of the threebody problem simulation.