Orbital physics for games

in Game Programming

Space, the final frontier and a recurring theme in games. The funny thing about space is that things rarely move in a straight line – then again, that all depends on your reference point. Gravity – and in particular orbital physics – forms an important aspect when it comes to movement in space. Many space games include some form of orbital mechanics or even depend on them as gameplay mechanic.

In this post I will talk how orbits in games can be implemented. In particular we will tackle some non-trivial orbits. This post will be very goal-focussed, and we will only discuss the underlying mathematics minimally, since it is primarily aimed to provide an easy reference to implement orbits yourself.

Additional sources

Before I start, I would like to point out two sources you can use alongside this post if you want more background information.

This first source is this GitHub repository which is an example project in Unity. This project was implemented by myself while figuring out the formulas myself. The code is MIT licensed, so you can copy the relevant components (Body and Orbit) into your own project if you want.

The second source is this Youtube channel. This channel was helpful in figuring out the last few things I couldn’t figure out myself and explains a few things more detailed than I will do in this post, so you are encouraged to watch the videos on the channel to supplement the information from this post.

What is an orbit?

(You can skip this section if you already understand the basic physics behind orbits.)

Orbits can be understood the best by using the thought experiment of Newton’s cannonball. If you drop a cannonball it will fall to the ground in a straight line because of gravity. If you fire a cannonball with some horizontal velocity, the cannonball will still fall to the ground, but in a trajectory depending on its initial velocity. If the initial velocity of the cannonball grows, the cannonball will land further and further away from its firing point. Because the Earth is round, the ground slowly curves away from the cannonball. If you fire the cannonball fast enough, the ground will curve away faster than the cannonball can fall, so the cannonball ends up in a state of endless falling towards the Earth in a circular, elliptic or hyperbolic path.

Newton's Cannon

The problem with full simulation

Games have the great property that they often have an update-loop. The game state is updated with relatively small time steps. This means that games are pretty good at doing physics simulations.

In particular, if we know for an object its current position \(\vec{p}\) and velocity \(\vec{v}\), we can calculate the position in a new frame as \(\vec{p} + \vec{v} \cdot \Delta t\). If we also know its acceleration \(\vec{a}\), we can calculate the new velocity as \(\vec{v} + \vec{a} \cdot \Delta t\). To compile, our update loop becomes something like this:

void Update(float dt)
{
    this.velocity += this.acceleration * dt;
    this.position += this.velocity * dt;
}

In our two-body orbit problem the acceleration is known, because the only force applied to the orbiting object (which we will call satellite in the general meaning of the word) is the gravity of the body it orbits. This acceleration is often referred to as \(\vec{g}\) and can be calculated (assuming that the mass of the body is much larger than the mass of the satellite) as \(\vec{g} = -\frac{G \cdot M}{r^2} \vec{r}\), where \(G\) is the gravitational constant, \(M\) the mass of the body, \(r\) the distance between the body and satellite, and \(\vec{r}\) the unit vector pointing from the body to the satellite. In code:

var r = this.position - this.body.position;
var rMag = r.magnitude; // or r.length, depending on your framework
var rNorm = r.normalised;
this.acceleration = -rNorm * (Constants.G * this.body.mass) / (rMag * rMag);

So with six lines of code, we have an update loop that is capable of simulating orbits about a body. End of story? Not really. Stepwise simulation – or numerical integration, which is what we are doing here – has a big problem: it’s unstable.

The particular integration method we have been using is called explicit Euler. The following image (shamelessly stolen from Wikipedia) shows this problem:

Euler method

At the initial state we use the current position and velocity to move forward. How much we move depends on the amount of time that has passed since the last frame. Even at a high fps, this will always be a discrete amount of time, and thus we will introduce a small error by moving in a straight line in the direction of the current velocity, since in a continuous simulation the velocity would be changing over time. In the second frame we therefore have a slightly perturbed position and velocity. While this error seems insignificant at first, it builds up over time, and with orbits in particular the satellite will slowly spiral away from the body.

There are different integration techniques that reduce this problem. Talking about these techniques is enough material for a post on its own. In this post we will look at a different stable solution. If you are interesting in starting with different integration techniques though, this article explains explicit Euler integration in more detail and also introduces RK4 integration, which is still not perfect, but better than explicit Euler and fairly easy to implement.

The beginnings: the circular orbit

The method we are going to discuss today is heavily based on Kepler’s Laws of Planetary Motion. Instead of using position and velocity to describe an orbit and the current progress along it, we will use Kepler’s elements.

Where position and velocity and intuitive and very likely to already be present in games, the Kepler elements are most likely to not be. Therefore we will start with very simple orbits and go from there.

Circular orbit

Above we see the most simple form of the two-body orbit problem: a satellite is rotating around a body at a constant distance \(r\). Independent of where the satellite currently is, we can easily move it to the next position by creating a rotation matrix that rotates the satellite around the body with the right amount. This immediately raises a question: what is the right amount?

Luckily the equations for completely circular orbits are very simple. From Newton’s Law of Universal Gravitation the following formula can be deduced: \[|\vec{v}| = \sqrt{\frac{G \cdot M}{r}}.\] Here \(G\) is the gravitational constant and \(M\) the mass of the body (which is assumed much larger than the mass of the satellite).

From the formula we can easily extract an angular velocity \(\omega\) (the angle in radians we rotate around the body) as well:
\[|\omega| = \frac{\sqrt{\frac{G \cdot M}{r}}}{2 \pi r}\].

// translational velocity
var v = Mathf.Sqrt(Constants.G * body.mass / r);
// translation this frame
var dp = v * Time.deltaTime;
// rotation this frame
var da = dp / (2 * Constants.PI * r);
// use a matrix transformation to rotate the satellite
this.RotateAround(parent.position, Vector3.up, da);

If we multiply this number by \(\Delta t\), we have our rotation. Of course we could also keep track of our current direction, update that each frame, and calculate or position from scratch.

Remark that since \(G\) is a constant and only influences the speed at which satellites orbit, we can freely choose it in our game. It is possible to remove the constant from the game entirely and balance the masses of the bodies only. It is also possible to disregard the exact speed formulas entirely and define the desired orbital period instead. To keep the velocities consistent, just keep in mind that the speed is halved if the radius is quadrupled.

Adding another dimension: inclination

Just the radius won’t cut it though. Given a body and a distance from that body, there are still infinitely many possible orbits that fulfil these parameters: we can rotate the orbital plane (the plane in which the orbit lies). That is: we can pick up one point on the orbit and move it up or down keeping it at the same distance from the body.

To fully define the circular orbit, we need to know the angle \(i\) the orbital plane makes with the plane that goes through the equator of the body. An inclination of 0 means we have an equatorial orbit, an inclination of 90 means we have a polar orbit, since the orbit passes over both poles, and an inclination of 180 means our satellite is on an equatorial orbit, but moves in opposite direction.

Still this is not enough information. This is why we also define the ascending and descending nodes. If we consider an orbit with a non-zero inclination, the orbit will intersect with the equatorial plane in two locations; the intersection where the satellite moves upwards (from the southern to the northern hemisphere) is called the ascending node, the opposite point the descending node.

Inclination

The ascending node will always lie on the equatorial plane and the distance to the body is already known as well. The only thing we need to define the position of the ascending node is its longitude \(\Omega\) (for bodies these are defined similarly to longitude on Earth).

Adapting our code to work with an inclined orbit is not difficult. The only thing we have to change is the vector we rotate around. We used Vector3.up before, but now we need a vector that is orthogonal to the orbital plane (note: the positive y-direction is assumed up and the equinox is assumed to point in the positive x-direction, analogue to the default Unity axes):

// we assume the inclination and longitude of ascending node to be given in radians
var axis = new Vector3(
    Mathf.Sin(inclination) * Mathf.Cos(longitudeOfAscNode),
    Mathf.Cos(inclination),
    Mathf.Sin(inclination) * Mathf.Sin(longitudeOfAscNode));

These are a bit eccentric: the elliptic orbit

So far so good, but we don’t live in a world where orbits are perfectly circular. In games, using circular orbits might be sufficient, either because the game mechanics don’t need more complicated orbits or circular orbits are sufficient to approximate most orbits. We will venture on to elliptic orbits however.

Remind that for a circle in two dimensions the points are given by \[\begin{pmatrix}\cos(\theta) \\ \sin(\theta)\end{pmatrix}.\]

For an ellipse we have a very similar structure. We will assume that the ellipse is aligned such that the long side is aligned with the x-axis. Further let \(a\) be the length of the semi-major axis (the distance between the centre and the furthest point on the ellipse) and \(b\) the length of the semi-minor axis (the distance between the centre and the closest point on the ellipse). The points are then given by \[\begin{pmatrix}a \cos(\theta) \\ b \sin(\theta)\end{pmatrix}.\]

Ellipse axes

With an elliptic orbit it is not sufficient to only know radius. This is why instead we generally speak of the height of the closest and furthest point in the orbit. The point in which the orbit is closest to the body is called the periapsis (for Earth orbits perigee is used and for Solar orbits perihelion; these three terms are often used interchangeably). The point directly opposite of that is the highest point in the orbit and is referred to as apoapsis (or apogee/aphelion).

Apoapsis/Periapsis

We will call the heights of the periapsis and apoapsis \(r_{per}\) and \(r_{apo}\) respectively. Because in an elliptic orbit the body will be one of the foci of the ellipse, we can use these numbers to calculate some important properties of the ellipse: \[a = \frac{r_{apo} + r_{per}}{2};\] \[e = \frac{r_{apo} – r_{per}}{r_{apo} + r_{per}};\] \[b = a \cdot \sqrt{1 – e^2}.\]

The original Kepler elements do not store the periapsis and apoapsis, but rather the periapsis and the eccentricity \(e\) of the orbit. In practice the periapsis and apoapsis are more intuitive, so we will assume those as given variables.

Just as with the ascending node, we need to know the location of the periapsis with respect to the body’s equinox. This is why the final parameter we need is the argument of periapsis. The argument of periapsis \(\omega\) (not to be confused with the angular velocity; sorry, I didn’t choose the symbols…) is the angle between the ascending node and the periapsis.

Argument of Periapsis

We can now easily calculate positions along the orbit, but there is one problem we still have to tackle. Where in circular orbits our speed is constant, this is not the case in elliptic orbits:

Elliptic orbit speed

The true anomaly \(\nu\) – the angle between the periapsis and the satellite’s current position as measured on the body – therefore does not change at a constant rate. Luckily there is a parameter that does change at a constant rate: the mean anomaly \(M\) (again: not to be confused with the mass). To convert between the true anomaly and mean anomaly, we need another anomaly to help us: the eccentric anomaly \(E\): the angle between the periapses and the satellite’s current position as measured in the centre of the ellipse. The eccentric anomaly in our case is also what we are interested in, because given the ellipse and the eccentric anomaly, we can calculate the satellite’s position.

Anomalies

The procedure below calculates the mean anomaly given the true anomaly:

float TrueToMeanAnomaly(float v)
{
    var e = eccentricity;

    var sinE = Mathf.Sin(v) * Mathf.Sqrt(1 - e * e) / (1 + e * Mathf.Cos(v));
    var cosE = (e + Mathf.Cos(v)) / (1 + e * Mathf.Cos(v));
    var E = Mathf.Atan2(sinE, cosE);

    return E - e * Mathf.Sin(E);
}

Remark that the intermediate variable E represents the eccentric anomaly.

The calculation in the other direction is more complicated. The reason is that we need to solve an equation that is not analytically solvable in the general case. This is why we use a numeric approximation procedure.

float MeanToTrueAnomaly(float M)
{
    var e = eccentricity;

    // make sure the mean anomaly lies in [-pi,pi]
    M %= 2 * Mathf.PI;
    if (M < -Mathf.PI)
        M += 2 * Mathf.PI;
    else if (M > Mathf.PI)
        M -= 2 * Mathf.PI;

    // calculate eccentric anomaly
    float E;
    if (M < 0)
        E = M - e;
    else
        E = M + e;

    var Enew = E;
    do
    {
        E = Enew;
        Enew = E + (M - E + e * Mathf.Sin(E)) / (1 - e * Mathf.Cos(E));
    } while ((Enew - E) > Mathf.Epsilon);
    E = Enew;

    // calculate true anomaly
    var sinv = Mathf.Sin(E) * Mathf.Sqrt(1 - e * e) / (1 - e * Mathf.Cos(E));
    var cosv = (Mathf.Cos(E) - e)/(1 - e * Mathf.Cos(E));
    return Mathf.Atan2(sinv, cosv);
}

The latter is the more interesting procedure. If we know our current anomaly, we can move it forward and calculate the new eccentric anomaly. Then we can calculate the new position of the satellite as well.

The only problem left is to decide how much to move the mean anomaly forward. Luckily there is a formula for that, just as with the velocity. The speed at which the mean anomaly changes is called the mean motion \(n\) and is calculated as \[n = \sqrt{\frac{G \cdot M}{a^3}}.\]

With all these elements, we can now write the update-loop for a satellite:

void Update(float dt)
{
    this.meanAnomaly += this.meanMotion * dt;
    var eccAnomaly = MeanToEccAnomaly(this.meanAnomaly);
    this.position = this.center + aVector * Mathf.Cos(E) + bVector * Mathf.Sin(E);
}

In this update loop, the MeanToEccAnomaly method is identical to the MeanToTrueAnomaly method, but skips the final lines in which the true anomaly is calculated and returns E instead. aVector and bVector are directional variants of the \(a\) and \(b\) parameters, i.e. the vector pointing in the direction of the semi-major and semi-minor axis with lengths \(a\) and \(b\) respectively.

Conclusion

In this admittedly longer than expected post I explained how to use Kepler elements to define orbits in games. I hope that especially the included code fragments will help using the explained mathematics in your own game.

If you spot any errors or have questions/feedback, feel free to leave a comment below.

Place comment

Your email address will not be published. Required fields are marked *