Finishing Our Gravity Simulator and “Solving” the Three Body Problem

Refreshing our memory

We are now ready to solve the Gravity Simulator Problem, reproduced below.

When we started planning a solution to this problem in a previous lesson, we wrote the following SimulateGravity() function to solve this problem.

SimulateGravity(initialUniverse, numGens, time)
    timePoints ← array of numGens+1 Universe objects
    timePoints[0] ← initialUniverse
    for every integer i between 1 and numGens
        timePoints[i] ← UpdateUniverse(timePoints[i-1], time)
    return timePoints

The preceding lesson’s primer on Newtonian dynamics led us to write UpdateUniverse() as follows.

UpdateUniverse(currentUniverse, time)
    newUniverse ← CopyUniverse(currentUniverse)
    for every body b in newUniverse
        oldA ← b.acceleration
        oldV ← b.velocity
        b.acceleration ← UpdateAcceleration(currentUniverse, b)
        b.velocity ← UpdateVelocity(b, oldA, time)
        b.position ← UpdatePosition(b, oldA, oldV, time)
    return newUniverse

To complete our gravity simulator, we will need to implement each subroutine called by UpdateUniverse(). We leave CopyUniverse() to the code alongs, since all this function does is produce a new Universe object by duplicating all the fields of an existing Universe object.

Updating the acceleration of a body

We will start with UpdateAcceleration(), which takes as input the Universe object currentUniverse as well as a Body object b, and returns an OrderedPair corresponding to the components of a vector representing the acceleration of b in the next time step.

This function passes its inputs to a subroutine ComputeNetForce() that returns the components of the net gravitational force acting upon b by all other bodies in currentUniverse. It then divides each of these components by the mass of b, according to the formulas below from our physics primer that rely upon Newton’s Second Law of Motion.

\begin{align*}
a_x(n+1) &= f_x(n+1) / m \\
a_y(n+1) &= f_y(n+1) / m
\end{align*}

The implementation of UpdateAcceleration() is shown below.

UpdateAcceleration(currentUniverse, b)
    accel ← OrderedPair object (0, 0)
    force ← ComputeNetForce(currentUniverse, b)
    accel.x ← force.x / b.mass 
    accel.y ← force.y / b.mass 
    return accel

Computing the net gravitational force acting on a body

To implement ComputeNetForce(), we range over all bodies in currentUniverse other than b, determine the components of the force exerted on b as a vector, and then add each of these components to the components of the net force. We rely on a subroutine ComputeForce(), which takes as input two Body objects b and b2 as well as the gravitational constant of currentUniverse and returns an OrderedPair corresponding to the components of the gravitational force of b2 acting on b, using this gravitational constant.

ComputeNetForce(currentUniverse, b)
    netForce ← OrderedPair object (0, 0)
    for every body b2 in currentUniverse other than b
        force ← ComputeForce(b, b2, currentUniverse.gravitationalConstant)
        netForce.x ← netForce.x + force.x
        netForce.y ← netForce.y + force.y
    return netForce

As for ComputeForce(), we start our function by determining the magnitude F of the gravitational force, which we know from our primer is given by the following equation, where the bodies have masses m₁ and m₂ and are separated by distance d, and where G is the gravitational constant.

F=\dfrac{G\cdot m_1 \cdot m_2}{d^2}

We can begin our function by establishing the initial x- and y-values of the force vector as equal to zero and then computing F, the magnitude of the gravitational force between the two input bodies.

ComputeForce(b, b2, G)
    force ← OrderedPair object
    force.x ← 0
    force.y ← 0
    d ← Distance(b.position, b2.position)
    F ← G * b.mass * b2.mass / (d^2)
    // to fill in

This function relies on a subroutine Distance(), which takes as input two OrderedPair objects corresponding to the positions of two objects and returns the distance between these objects’ positions.

Distance(p1, p2)
    deltaX ← p2.x – p1.x
    deltaY ← p2.y – p1.y
    d ← sqrt(deltaX^2 + deltaY^2)
    return d

To complete ComputeForce(), we need to allocate the gravitational force into components f_x and f_y given the positions of the bodies and the magnitude |F| of the force F between these bodies. As illustrated in the figure below, recall from the physics primer that if we know the angle (θ) indicating the direction from b to b2, then the following formulas hold:

\begin{align*}
|f_x| &= |F| \cdot \cos{\theta} \\
|f_y| &= |F| \cdot \sin{\theta} 
\end{align*}

Using this idea, we can determine the angle formed between b and b2 using a subroutine Angle(), and then apply the above equation to determine the components of the force vector, as shown by the following completed ComputeForce() function.

ComputeForce(b, b2, G)
    force ← OrderedPair object
    force.x ← 0
    force.y ← 0
    d ← Distance(b.position, b2.position)
    F ← G * b.mass * b2.mass / (d^2)
    theta ← Angle(b, b2)
    force.x ← F * Cosine(theta)
    force.y ← F * Sine(theta)
    return force

Although most programming languages have built-in functions for sine and cosine, writing the Angle() subroutine will require some involved trigonometry. Yet we can determine the cosine and sine of θ without ever having to compute θ. To do so, note that we have already calculated the distance d between the two bodies. We also know that d is the length of the hypotenuse of a right triangle whose base has length equal to the positive difference between the x-coordinates of the bodies’ positions, which we denote deltaX. As a result, cos(θ) is equal to the ratio deltaX/d, and a similar argument gives us that sin(θ) is equal to the ratio deltaY/d.

We therefore will update ComputeForce() below to first determine deltaX and deltaY. Then, we will replace Cosine(theta) and Sine(theta) with deltaX/d and deltaY/d, respectively. This approach removes the need for an Angle() subroutine altogether.

ComputeForce(b, b2, G)
    force ← OrderedPair object
    force.x ← 0
    force.y ← 0
    d ← Distance(b.position, b2.position)
    F ← G * b.mass * b2.mass / (d^2)
    deltaX ← b2.position.x – b.position.x
    deltaY ← b2.position.y – b.position.y
    force.x ← F * deltaX / d
    force.y ← F * deltaY / d
    return force

Updating a body’s velocity and position

Now that we have written UpdateAcceleration(), we will write UpdateVelocity() and UpdatePosition(). Fortunately, we have equations for updating the components of velocity and position from our physics primer.

\begin{align*}
v_x(n+1) &= \tfrac{1}{2}\,\big(a_x(n) + a_x(n+1)\big)\,t + v_x(n) \\
v_y(n+1) &= \tfrac{1}{2}\,\big(a_y(n) + a_y(n+1)\big)\,t + v_y(n)
\end{align*}

\begin{align*}
p_x(n+1) &= \tfrac{1}{2}\,a_x(n)\,t^2 + v_x(n)\,t + p_x(n) \\
p_y(n+1) &= \tfrac{1}{2}\,a_y(n)\,t^2 + v_y(n)\,t + p_y(n)
\end{align*}

We now simply need to translate these equations into pseudocode, as shown below.

UpdateVelocity(b, oldA, t)
    newA ← b.acceleration   // must already be a(n+1)
    newV ← OrderedPair object
    newV.x ← b.velocity.x + 0.5 * (oldA.x + newA.x) * t
    newV.y ← b.velocity.y + 0.5 * (oldA.y + newA.y) * t
    return newV

UpdatePosition(b, oldA, oldV, t)
    newP ← OrderedPair object
    newP.x ← b.position.x + oldV.x * t + 0.5 * oldA.x * t^2
    newP.y ← b.position.y + oldV.y * t + 0.5 * oldA.y * t^2
    return newP

It’s time to code

Our gravity simulator is now complete, and all that remains is to implement it. Please join us in the following code alongs to learn how to implement object-oriented programming and generate some beautiful gravity simulations that will save the residents of Trisolaris once and for all and knock Newton’s stockings off.