Building a Gravity Physics Engine in Go

Learning objectives

In this code along, we will build a physics engine to simulate the force of gravity and how it affects the motions of a collection of heavenly bodies, solving the problem below whose solution we discussed in a language-neutral fashion in the core text. We will then run this simulation on the moons of Jupiter as well as a collection of three-body simulations that exhibit beautiful patterns.

All of our physics engine functions will live in gravity.go. We will implement them from the top down, starting with SimulateGravity() and working our way down to the individual force and kinematic calculations.

Code along summary

Copying a universe

Before writing the simulation, we need to address a subtle but important issue: deep copying. At each time step, we want to create a new Universe that is an independent copy of the previous one. We cannot simply write newUniverse = currentUniverse, because Universe contains a slice (bodies), and in Go, assigning a slice only copies the slice header — both variables would still point to the same underlying array. Modifying one would silently corrupt the other.

Instead, we write explicit copy functions. CopyBody() copies all the fields of a single Body into a new one. Since Body contains only scalar values and OrderedPair structs (which are themselves plain value types with no slices or pointers), a field-by-field copy is sufficient.

// CopyBody returns a new Body whose fields are identical to b.
func CopyBody(b Body) Body {
	var b2 Body
	b2.name = b.name
	b2.mass = b.mass
	b2.radius = b.radius
	b2.red = b.red
	b2.green = b.green
	b2.blue = b.blue
	b2.position.x = b.position.x
	b2.position.y = b.position.y
	b2.velocity.x = b.velocity.x
	b2.velocity.y = b.velocity.y
	b2.acceleration.x = b.acceleration.x
	b2.acceleration.y = b.acceleration.y
	return b2
}

CopyUniverse() copies the scalar fields of a Universe directly, then creates a fresh slice of bodies and copies each body individually using CopyBody().

// CopyUniverse returns a deep copy of currentUniverse.
func CopyUniverse(currentUniverse Universe) Universe {
	var newUniverse Universe
	newUniverse.gravitationalConstant = currentUniverse.gravitationalConstant
	newUniverse.width = currentUniverse.width

	// Make a brand-new slice; never assign slices directly
	numBodies := len(currentUniverse.bodies)
	newUniverse.bodies = make([]Body, numBodies)
	for i := range newUniverse.bodies {
		newUniverse.bodies[i] = CopyBody(currentUniverse.bodies[i])
	}

	return newUniverse
}

SimulateGravity

SimulateGravity() is our top-level function. It allocates a slice of numGens + 1 universes, sets the first element to initialUniverse, and fills in each subsequent element by calling UpdateUniverse().

// SimulateGravity returns a slice of numGens+1 Universe snapshots,
// starting from initialUniverse, with each step separated by time seconds.
func SimulateGravity(initialUniverse Universe, numGens int, time float64) []Universe {
	timePoints := make([]Universe, numGens+1)
	timePoints[0] = initialUniverse
	for i := 1; i < numGens+1; i++ {
		timePoints[i] = UpdateUniverse(timePoints[i-1], time)
	}
	return timePoints
}

UpdateUniverse

UpdateUniverse() computes one time step. It first deep-copies the current universe so that we can safely read from the original while writing the new values. For each body, it updates the acceleration, velocity, and position in order, taking care to save the old acceleration and old velocity before they are overwritten, since the update equations for velocity and position both depend on those old values.

// UpdateUniverse returns a new Universe representing one time step forward.
func UpdateUniverse(currentUniverse Universe, time float64) Universe {
	newUniverse := CopyUniverse(currentUniverse)

	for i, b := range newUniverse.bodies {
		oldAcceleration, oldVelocity := b.acceleration, b.velocity
		newUniverse.bodies[i].acceleration = UpdateAcceleration(currentUniverse, b)
		newUniverse.bodies[i].velocity = UpdateVelocity(newUniverse.bodies[i], oldAcceleration, time)
		newUniverse.bodies[i].position = UpdatePosition(newUniverse.bodies[i], oldAcceleration, oldVelocity, time)
	}

	return newUniverse
}

Notice that we iterate over newUniverse.bodies but read gravitational forces from currentUniverse. This ensures all bodies’ forces are calculated from the same snapshot of positions, rather than using a mix of old and already-updated positions.

Computing acceleration

The acceleration of each body is determined by the net gravitational force acting on it from all other bodies. By Newton’s second law, F = ma, so a = F/m. UpdateAcceleration() computes the net force vector and divides it componentwise by the body’s mass.

// UpdateAcceleration returns the acceleration of body b due to all other bodies
// in currentUniverse, computed via Newton's second law F = ma.
func UpdateAcceleration(currentUniverse Universe, b Body) OrderedPair {
	var accel OrderedPair
	force := ComputeNetForce(currentUniverse, b)
	accel.x = force.x / b.mass
	accel.y = force.y / b.mass
	return accel
}

ComputeNetForce() loops over all bodies in the universe and sums the individual gravitational force contributions, skipping the body b itself.

// ComputeNetForce returns the net gravitational force vector acting on b
// from all other bodies in currentUniverse.
func ComputeNetForce(currentUniverse Universe, b Body) OrderedPair {
	var netForce OrderedPair // starts at (0, 0)

	for i := range currentUniverse.bodies {
		if currentUniverse.bodies[i] != b {
			G := currentUniverse.gravitationalConstant
			currentForce := ComputeForce(b, currentUniverse.bodies[i], G)
			netForce.x += currentForce.x
			netForce.y += currentForce.y
		}
	}

	return netForce
}

The force between two bodies is given by Newton’s law of universal gravitation: F = G m₁ m₂ / d², where d is the distance between them. The force acts along the line connecting the two bodies, so we break it into x and y components using the direction vector (Δx/d, Δy/d).

// ComputeForce returns the gravitational force vector acting on b due to b2,
// using G as the gravitational constant.
func ComputeForce(b, b2 Body, G float64) OrderedPair {
	var force OrderedPair

	d := Distance(b.position, b2.position)
	if d == 0.0 {
		return force // avoid division by zero
	}

	F := G * b.mass * b2.mass / (d * d)

	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
}

// Distance returns the Euclidean distance between two positions.
func Distance(p1, p2 OrderedPair) float64 {
	deltaX := p1.x - p2.x
	deltaY := p1.y - p2.y
	return math.Sqrt(deltaX*deltaX + deltaY*deltaY)
}

Updating velocity

We update the velocity of a body using the average acceleration method, which averages the acceleration at the beginning and end of the time step. This gives a more accurate numerical integration than simply using the acceleration at one endpoint. The update formula is:

v(t + Δt) = v(t) + 0.5 * (a(t) + a(t + Δt)) * Δt

// UpdateVelocity returns the updated velocity of body b after one time step.
// oldAcceleration is b's acceleration at the start of the step; b.acceleration
// has already been updated to the end of the step.
func UpdateVelocity(b Body, oldAcceleration OrderedPair, time float64) OrderedPair {
	var currentVelocity OrderedPair
	currentVelocity.x = b.velocity.x + 0.5*(b.acceleration.x+oldAcceleration.x)*time
	currentVelocity.y = b.velocity.y + 0.5*(b.acceleration.y+oldAcceleration.y)*time
	return currentVelocity
}

Updating position

We update the position using constant-acceleration kinematics, using the body’s velocity and acceleration at the start of the time step:

p(t + Δt) = p(t) + v(t) * Δt + 0.5 * a(t) * Δt²

// UpdatePosition returns the updated position of body b after one time step.
// oldAcceleration and oldVelocity are b's values at the start of the step.
func UpdatePosition(b Body, oldAcceleration, oldVelocity OrderedPair, time float64) OrderedPair {
	var pos OrderedPair
	pos.x = b.position.x + oldVelocity.x*time + 0.5*oldAcceleration.x*time*time
	pos.y = b.position.y + oldVelocity.y*time + 0.5*oldAcceleration.y*time*time
	return pos
}

Running the simulation

With all the physics functions in place, we can run our simulation. From inside your go/src/gravity directory, execute the following command to simulate the moons of Jupiter for 1000 time steps of 60 seconds each, rendered on a 1500 × 1500 canvas and saving every 10th frame.

go run *.go jupiterMoons 1000 60 1500 10

You can also run three-body simulations that exhibit remarkable periodic orbits. Here are some examples to try.

go run *.go butterfly 100000 0.0001 1500 300
go run *.go dragonfly 100000 0.0001 1500 300
go run *.go figureEight 100000 0.0001 1500 300
go run *.go goggles 100000 0.0001 1500 300
go run *.go moth 100000 0.0001 1500 300
go run *.go yarn 100000 0.0001 1500 300

Each command produces an animated GIF in the output/ directory. The arguments are, in order: the scenario name, the number of time steps, the time interval per step (in seconds), the canvas width in pixels, and the drawing frequency (one frame is saved every n time steps).

Congratulations! You have now implemented a complete gravity simulator in Go. The combination of Newton’s law of universal gravitation with the average acceleration and kinematic position update equations is powerful enough to accurately model the orbits of real celestial bodies as well as the beautiful choreographed three-body solutions discovered by mathematicians and physicists.