The Keys to the Universe

Okay, I lied a little bit with the title. But I have been playing around with universe simulation stuff (or, to use the formal title, N-body gravitational simulation) recently and wanted to share what I’ve learned with the void. This is going to be a long and technical post, so be warned!

One of the ideas I had for the latest GMTK game jam was to make a countdown until asteroid impact game. I quickly ruled it out because of the complicated physics calculations that would need, and instead focused on something else. It turns out that was the right call, having just spent a few days building what I have done.

But, I was still curious about the concept, so I took it upon myself to explore this region for a few days. It turns out that I’m not alone, which is probably not such a bad thing, because I wouldn’t have made it very far otherwise. I think I’ve spent more time searching, watching videos and reading articles than developing solutions over the last handful of days.

In the hopes of saving myself future headaches, and anyone else that stumbles upon this corner of the web, I’m writing down the lessons I’ve learned so far.

But first, a satisfying, albeit messy looking results video of these efforts so far

I mean, there were supposed to be two extra moons orbiting those planets.

Getting Things Moving

How do we get objects moving about under their own force? For that, we need a few basic pieces of information about each object. For my system, I’ve called these objects a Body – whether it’s a star, planet, moon, asteroid, they all follow the same basic model here (I do love to complicate matters.)

To start with, all we need is a position and a mass. Well, velocity too, but I have that calculated during the initialisation stages.

Euler vs Verlet Integration

You may have heard about explicit Euler, semi-implicit Euler, position Verlet, velocity Verlet, or even Runge-Kutta (RK4) methods for moving objects in a game. Even if you haven’t, you’re almost certainly familiar with the following:

C#
void UpdatePosition(float deltaTime)
{
    Acceleration = CalculateAcceleration();
    Position += Velocity * deltaTime;
    Velocity += Acceleration * deltaTime;
}

That is an example of the explicit euler method. It’s also a terrible approach for simulating physics in this scenario, because of compounding approximation errors. Using it for a character controller might be okay, but when calculating orbital movements (which can involve hundreds of millions of the same calculation), the errors stack up quickly. This is where the concept of numerical stability comes into play – a term that confused me when I first heard it (how can numbers be unstable?!)

The approaches I was considering were:

  1. Explicit Euler
  2. Semi-Implicit Euler
  3. Position/Velocity Verlet
  4. Runge Kutta

It’s worth noting that there are benefits and costs to each approach, it’s not as simple as worst to best. Explicit Euler works great for many common game development needs and the extra computation isn’t worthwhile. This is just one of those situations where it doesn’t offer the right degree of accuracy, so a more expensive computation is required. RK4 does things I don’t even understand, so don’t quote me on any of that.

Given that RK4 computes it’s data 4 times to get it’s accuracy, it’s an expensive option to pick, so I went with the next best thing. I tried position verlet first, but settled on velocity verlet as it was easier to grasp (and makes the initialisation section described later a bit easer)

The concept is straightforward enough (see the actual formulas on the Velocity Verlet Wikipedia article, or my code implementation below):

  1. Calculate the new position using the current velocity and acceleration
  2. Calculate the new acceleration using your selected method (forces, user input, etc)
  3. Calculate the new velocity using the current velocity, current acceleration and new acceleration

For something that seems so simple, it is remarkably stable. Once I had it up and running, my error drift (the unintentional value change between start and end) over a few thousand orbits, or a several hundred million calculations, was around ~10-11 (which I still can’t quite believe seeing) – don’t ask me how it is stable. I spent half an hour getting completely confused over phase-space areas and 6-dimension something that explains it. I just know that it works.

My implementation for this is about as simple as it can get:

C#
public void Tick(SimulationState state)
{
	foreach (var body in state.Bodies)
	{
		body.NewPosition = body.Position + body.Velocity * SimulationTimestep + (0.5 * SimulationTimestep * SimulationTimestep * body.Acceleration);

		body.NewAcceleration = SimulationPhysicsHelper.CalculateAccleration(state, body, body.NewPosition);

		body.NewVelocity = body.Velocity + 0.5 * SimulationTimestep * (body.Acceleration + body.NewAcceleration);
	}
}

Calculating Gravity

We can thank Newton for this one. And probably many other people too. To calculate the forces applied, we need to use Newton’s law of universal gravitation, which gives us the following formula:

F=GMmr2F = G\frac{Mm}{r^2}

(where G is the gravitational constant – which you could just make up for a game, M and m are the respective body masses, and r is the distance between the two bodies)

Throw in a quick F=ma and juggle the formulas a little to get the acceleration, which is that mystical third property we’re calculating. Of course, F here is just a number – to translate that into a vector, we simply multiply the direction towards the other object.

I loop over every object in the simulation, add up all of the force vectors, then divide that by the mass to get a total acceleration value

C#
public static Vector2D CalculateAccleration(SimulationState state, Body body, Vector2D position)
	=> CalculateForces(state, body, position) / body.Mass;

public static Vector2D CalculateForces(SimulationState state, Body body, Vector2D position)
{
	Vector2D total = Vector2D.zero;

	foreach (var otherBody in state.Bodies)
		if (body != otherBody)
			total += CalculateForceFrom(otherBody, position, body.Mass);

	return total;
}
public static Vector2D CalculateForceFrom(Body body, Vector2D position, double mass)
{
	var direction = body.Position - position;
	var distanceSqr = direction.sqrMagnitude;

	var force = GravitationalConstant * (body.Mass * mass) / distanceSqr;

	return direction.normalized * force;
}

Initialisation

Great – I can accurately move items around. But they just fall into the sun. To get these planets actually orbiting, we need to give them an initial push. Gravity will then act as a sort-of centripetal force to maintain the orbit, no extra effort needed.

Ironically, this is one of the harder parts to figure out. The jump-start calculations that I needed once, rather than the (altogether simpler) update formulas that run millions of times.

That’s because I love complexity, potentially at the cost of my own sanity. Other solutions will often take a shortcut, and say “this moon orbits that planet, which orbits that star.” I didn’t want that. I wanted to place my objects in space, define it’s mass and hit play.

Thankfully, someone (or someones) very clever have once again provided the perfect formula for determining the circular orbit velocity:

v=G(M+m)av=\sqrt{\frac{G(M+m)}{a}}

(where G is the gravitational constant, M and m are the masses of the two bodies, and a is the length of the “semi-major axis” – or simply the radius for a perfect circle)

C#
public static Vector2D CalculateSphericalOrbitalVelocity(Body dominantBody, Body orbitingBody)
{
	var direction = orbitingBody.Position - dominantBody.Position;
	var tangent = new Vector2D(-direction.y, direction.x).normalized;

	var distance = direction.magnitude;
	double masses = dominantBody.Mass + orbitingBody.Mass;

	return tangent * Math.Sqrt(GravitationalConstant * masses * / distance);
}

Using this equation gave me over 4000 stable orbits in one of my tests (once I figured out the below issues, that is). Work out the tangent direction between the two bodies, and you’ve got your starting velocity! Now it’s just time to plug it in and go!

Orbital Parent Selection

Of course, if you’re making a hierarchy-less simulation, getting the right velocity means determining what each object should be orbiting. It’s worth noting that this is only for the sake of initialisation, to determine the appropriate starting velocity of a body. This concept of orbital parent doesn’t get saved with the body.

Simply using the most massive body works when you only have planets orbiting the sun, but stops working as soon as you want moons orbiting planets. Turned out, this took a couple of different approaches all merged into one to figure out. For this, I’ve mentally coined the term orbital parent – there’s no hierarchy, so the moon can leave whenever it wants (which it frequently did in my simulations once I added many more planets.)

For this, I relied on an equation known as the Hill Sphere radius to determine the radius of influence for bodies orbiting around a more massive object (such as planets orbiting the sun), giving an approximate radius in which the orbiting body can maintain its own satellites, despite the larger force applied by it’s orbital parent.

r=aM3(M+m)r=a\sqrt{\frac{M}{3(M+m)}}

(where M is the mass of the more massive body (i.e. sun) and m is the mass of the smaller body (i.e. planet), a is that pesky semi-major axis, or radius for a perfect circle)

While the equation labels itself as approximate (it works best for a two-body system rather than an N-body system), this gives surprisingly good results, but with one important caveat – I have to measure the candidate of the candidate body and it’s orbital parent. This means that for any given candidate, I need to find it’s orbital parent too – a recursive problem of sorts. It looks a bit like this:

  1. Loop through existing bodies to find object more massive than current body
  2. Determine orbital parent of existing candidate parent
  3. Calculate hill sphere radius for candidate and it’s parent
  4. If within radius, return candidate parent as current body’s orbital parent

There’s a problem with this approach – you can’t calculate the hill sphere radius of the Sun. While it’s without a doubt the most massive object, it doesn’t have an orbital parent in which its own gravity dominates (well, not in this simulation anyway, we’ll keep it strictly limited to one solar system please)

For that, I use a simple fallback mechanism: if the hill sphere method finds no suitable parent, we’ll go back to Newton’s equation to determine the force applied between all objects, and return the largest force. Not quite so elegant, but works pretty well.

C#
private Body GetInfluencingBody(SimulationState state, Body newBody)
{
	double largestForce = 0;
	Body largestForceBody = null;

	foreach (var body in state.Bodies)
	{
		if (body == newBody)
			continue;

		if (body.Mass <= newBody.Mass)
			continue;

		//Attempt to use hill sphere radius first
		//Hill sphere is between the body and _it's_ parent, not the new body (between planet and sun, not moon and planet)
		Body orbitalParent = GetInfluencingBody(state, body);

		if (orbitalParent != null)
		{
			double radius = PhysicsCalculations.CalculateHillSphereRadius(orbitalParent, body);

			double distance = Vector2D.Distance(body.Position, newBody.Position);

			if (distance < radius)
				return body;
		}

		//If that doesn't work, fall back to largest gravitational force
		var forceMagnitude = PhysicsCalculations.CalculateForceFrom(body, newBody.Position, newBody.Mass).magnitude;

		if (forceMagnitude > largestForce)
		{
			largestForce = forceMagnitude;
			largestForceBody = body;
		}

		//Todo: Return multiple candidates if they exist
	}

	return largestForceBody;
}

Excellent! We now have everything needed to put together more complex simulations.

Floating Point Precision

I haven’t heard this talked about, but I immediately ran into an issue that sank most of the weekend trying to resolve. Mainly because I wasn’t sure if the integration was incorrect, my scale was off, or if there was some other bug in the code. All I knew at the time was that I was getting a gradual increase of orbital distance in my two body test case. After only 100 orbits, my planet was noticeably further out of position.

Orbit 1: 99.97115 -> 100.2975
Orbit 101: 107.9131 -> 109.9531

I felt like I was going around in circles trying to fix the problem – I even started the whole project again, assuming I’d simply missed a bug. Thankfully the concept of floating point precision was high on my list. Since I was dealing with a bunch of large and small numbers (even scaled down by 1024 or so, the mass of the Sun is a pretty big number, and the gravitational constant I chose was also slightly scaled version of the real thing at 6×10-5 or so)

Unity works using floats, mainly. Anyone that tries to edit things according to a grid has likely seen their perfect rounded numbers change from 100 to 99.99999999. After trying a few things (including a complete re-write), the next step was to try using doubles. That involved creating my own implementation of Vector2, replacing the Mathf functions with Math, and a few direct replacements of values and types within the rest of the code.

A fair bit of work, that I wasn’t too keen on trying, but boy am I glad I did!

Orbit 10: 149.600006103496 -> 149.600006400084
Orbit 4800: 149.600006102133 -> 149.600006398786

That’s basically nothing!

The reason for this is quite straightforward – double offers 8 bytes for storing the value, whereas float is just 4, giving me from roughly 6-9 digits up to roughly 15-17 digits of precision. There’s probably a cost to doing this, but since my approach was simply broken without it, I had no choice.

My initial approach using floats involved multiplying a very small number with a very large number, and repeating that thousands of times per second. It’s no wonder that errors started creeping in. That said, I do wonder if the new simpler values that I have would work better, but I’m not in a hurry to test that theory at the moment.

Eccentric Orbits

Not all orbits are perfectly circular. In fact, I’d be surprised if any were. I definitely wanted to get some more, er, “oval” shaped orbits in there. This is where I came across Orbital Eccentricity. With one extra number and a tiny change to my earlier formula, I can add the ability to have planets orbit with eccentric orbits. The formula goes like:

v=G(M+m)(1+e)av=\sqrt{\frac{G(M+m)\color{red}{(1+e)}}{a}}

(where e is the eccentricity)

While a simple change, this formula does carry an important limitation which I’ve expanded on below – that is, it only calculates the velocity from the closest point of orbit, so the orbit won’t get closer than it already is. That’s fine for a next step in the simulation, though.

e is any non-negative value that describes what type of trajectory the body will have.

When e=0, the body has a circular orbit (as you can see in the above formula, multiplying by (1+0) makes no change to the value, so would be the same as the earlier circular orbit formula.

For e>0 && e<1, this is where the magic of the eccentric orbit comes into play. Higher values produce a more stretched oval, lower values tend back towards the circular orbit.

While e>1, the body is on a hyperbolic trajectory rather than an orbit – in short, the gravity of its parent affects its direction but it’s going fast enough to escape the gravitational pull of the parent.

And finally, e=1 is what they call a parabolic trajectory. It’s a special case, where the body approaches infinity, it’s velocity approaches zero so it may never technically escape the gravitational pull of it’s parent. For all practical purposes in any game we’d make, it is as good as the above hyperbolic trajectory though.

A simple modification to the circular orbit to allow this to function:

C#
public static Vector2D CalculateEccentricOrbitalVelocity(Body dominantBody, Body orbitingBody, double eccentricity)
{
	if (eccentricity < 0)
		throw new ArgumentException("Value must not be negative", nameof(eccentricity));

	var direction = orbitingBody.Position - dominantBody.Position;
	var tangent = new Vector2D(-direction.y, direction.x).normalized;

	var distance = direction.magnitude;
	double masses = dominantBody.Mass + orbitingBody.Mass;

	return tangent * Math.Sqrt(GravitationalConstant * masses * (1 + eccentricity) / distance);
}

Now, doesn’t that make for a much more interesting simulation.

Where I’m At

So, I have a simulation that’s capable of handling circular and basic eccentric orbits, updating positions in a much-less-error-prone method using Newton’s law of universal gravitation and velocity verlet integration.

I can produce a stable two-body circular orbit with negligible numerical drift over thousands of orbits and hundreds of millions of step calculations.

I can simulate over a hundred bodies all interacting and influencing each other, calculating initial orbit velocities using a hierarchy-less, data-driven approach.

Not too bad for a few days work.

What’s Next?

That’s all I have for now. I have proven, to myself if nobody else, that simulating orbits without explicit hierarchy is doable. I have a few thoughts on future directions, but for now, the formulas involved are getting a bit too far above my head.

Eccentric Orbits

The eccentric orbit formula currently assumes the starting position as the minimum distance from the star. I’d like to be able to change the exact position on the elliptical orbit to allow for a bit more control over how it first works.

Binary Orbits

Putting two stars of equal mass next to each other just means they start with no velocity (because neither is more massive than the other). Equally, giving even the tiniest amount more mass means that one is launched with an orbit, rather than both stars orbiting each other. I’d like to add in support for calculating the barycentre of objects that are massive enough and close enough to influence each other, in binary-like orbits. The follow up from this (if you read the todo in my earlier code example, you’ll have seen this coming) would be to have orbits around binary orbits, whereas right now, it only supports finding a single orbital parent.

Pre-Warming

This one is pretty easy in concept, just run the simulation for a few orbits, but it may be computationally expensive. Still, having the orbits already in motion would make for a better view of the game when/if it unfolds further.

Trajectory Lines

Again, a simple sounding feature, but potentially quite expensive. I haven’t given it much thought though, other than “this would be a good feature to add”

Collisions

Naturally. Objects don’t just phase through each other (well, they do in my simulation) – it currently leads to interesting effects like asteroids getting ejected from the solar system. In reality, this would need to disrupt both orbits, possibly fragment the bodies if I really wanted to go that far, and a whole bunch of fun sounding physics things to go with it. If I want a game out of that, this is high on the list to look at next.

Performance

I didn’t talk too much about performance, but lets say it’s a O(n2) kind of problem. That is, for every extra body I add to the simulation, the computational cost goes up quadratically. With some 140-ish bodies in the simulation, my framerate was halved from 60 (v-sync capped) to 30

That’s because each of the 140 bodies were checking against the other 139 for force calculations. That’s a lot of sums in a short time period. I have some simple tricks to use, some algorithmic, and some to take advantage of technology to improve performance..

First, I can skip lower mass objects if they’re sufficiently smaller or far enough away that their force would be negligible. This might need a bit of trial and error in dialling the numbers, but the sun probably doesn’t need to care about asteroids pulling on it.

My formulas above are not exactly written in a performant manner (for example, I’m running 0.5 * SimulationTimestep * SimulationTimestep every iteration, and there are many, many more examples.) Simplifying many of these equations may make them less readable, but would reduce unnecessary calculations. For purposes of my experiment, I’ve left them in a longer form.

This is also a great use-case for multi-threading. Positions are calculated before any bodies are moved, so would benefit from multi-core processing that’s a lot more available in modern-day processors. I’ve heard about Unity ECS/burst and this may be an area to look into in the future.

I’ve also come across Barnes-Hut in my research, which appears to be a way to group objects into regions, so calculations use a single, averaged force and distance instead of working out each object individually – much more promising when dealing with many objects further away.