# Interactive Astrodynamics

Kepler's laws of planetary motion describe the elliptical path of a object in a stable orbit.
This includes the paths of objects like the Moon and the International Space Station around the Earth, or the path of the Earth around the Sun.

Newton used his own laws to derive Kepler's, showing that celestial objects moved according to the same physical rules as those closer to home.
Indeed, it is trivial to [implement an _n_-body simulation](https://youtu.be/Ouu3D_VHx9o?t=403) of gravity using Newton's law of gravitation.
However, this approach works only for very simple scenarios.
Physics calculations in games are done with floating-point arithmetic, of which a limitation is that [numbers may only be approximately represented](https://www.youtube.com/watch?v=PZRI1IfStY0).

Limited precision poses a subtle and hard to solve problem which is present in any iterative algorithm, such as an _n_-body simulation based on universal gravitation.
As part of its operation, the state at the end of one discrete [frame](https://gameprogrammingpatterns.com/game-loop.html) forms some or all of the input for the next.
This means that error due to the lack of precision accumulates over time, resulting in larger discrepancies from reality the longer the simulation is run for.

Unfortunately, floating points are a fact of life in game physics, and we are not going to escape them here.
So how does Kerbal Space Program do it?

![ksp.jpg](ksp.jpg)

Kepler's laws make a simplification known as the two-body problem.
More complex systems, such as Sun-Earth-Moon or the entire Solar System can be modelled as a composition of several two-body systems, rather than the _n_-body systems they really are.
This conveniently sidesteps the requirement to use approximate forces and therefore accumulate error over time.
The two body problem describes the motion of one body around another, or more precisely a common barycentre.
In the case of one massive and one tiny object, the barycentre is approximately in the same position as the centre of the larger object and can therefore be ignored.
In a two-body system, an orbit is perfectly stable and cyclical.
A body's position can thus be calculated directly as a function of time.

The ability to calculate positions directly as a function of time has some extremely valuable properties.
One of the many challenges faced in deep space exploration is the sheer duration required by such voyages - the Voyager probes were launched nearly half a century ago, and even reaching Mars takes the better part of a year.
To be of any real use, then, interactive astrodynamic simulations need a level of time control.
Kerbal Space Program has a time warp feature, allowing players to navigate missions taking years of simulated time in an evening of real time.

This would not be possible with an iterative physics simulation, as the accuracy of the result is directly dependent on the resolution of the time step.
Games have to produce several frames per second to both sell the illusion of motion and stay responsive, meaning each frame must be computed in a small fraction of a second.
Depending on the platform, target frame rates can vary from as low as 30Hz, for a mobile device preserving its battery, to 120Hz and beyond for competitive or virtual reality applications.
This means most games have about 10ms to compute each frame.

Processing limitations are of particular interest in games, which are severely constrained in the amount of processing they are able to complete in a single rendered frame.
In practice, each sub-system comprising a game has significantly less than 10ms, since it must use its budget conservatively to allow other sub-systems to use the same resources.
An astrodynamics physics simulation may be allowed only a single millisecond or less in which to calculate the paths of objects, due to constraints imposed by resourcing requirements of other sub-systems.

Simulating physics to a high degree of accuracy with an iterative model requires high-resolution time steps.
Interactive physics simulations will often run at many times the rendered framerate to maintain stability.
For example, simulating an [internal combustion engine](https://www.youtube.com/watch?v=RKT-sKtR970) requires extremely high simulation frequencies to be able to accurately produce sound.
To maintain accuracy at a higher rate of passage of time, the simulation frequency has to increase with it.
Kerbal Space Program's highest level of time warp is 100,000x normal speed, at which a physics simulation would require computation of 100,000x the time steps per rendered frame to remain stable.

If we instead calculate an object's position from a timestamp, the calculus changes completely.
The simulation is no longer constrained by the amount of simulated time the CPU can process per unit real-time.
Instead, we are able to jump directly to the time of our choosing, at whatever rate we like, and determine the position of the object afresh every time.
The object's position ceases to be dependent on its previous position, and the cost of our calculation becomes unbound from the rate of passage of simulated time, so long as we maintain a constant frame rate.
We can even run time backwards - from the simulation's perspective, it does not care a jot.

---

Planets move in elliptical paths, and [ellipses are surprisingly complex](https://www.chrisrackauckas.com/assets/Papers/ChrisRackauckas-The_Circumference_of_an_Ellipse.pdf).
Part of the determination of a planet's position as a function of time is the calculation of the [eccentric anomaly](https://en.wikipedia.org/wiki/Eccentric_anomaly) from the mean anomaly, which is related by Kepler's equation.

![keplers-equation.png](keplers-equation.png)

Having no closed-form solution for _E_, the equation must be solved with numerical methods.
A common choice is [Newton-Raphson](https://en.wikipedia.org/wiki/Newton%27s_method), a general-purpose root-finding algorithm which converges on successively better approximations through an arbitrary number of iterations.
It is straightforward to understand and implement, but it has some major flaws for real-time physics simulation.

For low-eccentricity orbits - those that are nearly circular - Newton-Raphson converges quickly in just a handful of iterations.
However, when presented with an eccentric orbit, such as those of the [Juno probe](https://science.nasa.gov/mission/juno/), NM takes exponentially more iterations to resolve to a high-accuracy solution.
With a lower number of iterations, the calculated value for _E_ is erratic and inaccurate, completely unsuitable for a stable simulation.

This presents a massive problem in a real-time system like a game.
To be fit for purpose a physics simulation needs to be stable in all reasonable cases the player might encounter, so a Newton-Raphson-based approach would need to be given the resources to run as many iterations as needed to determine a suitably accurate value for _E_.
Since the behaviour is of an exponential increase in computation as orbits become increasingly eccentric, the worst-case performance is many orders of magnitude worse than the best- or even average-case performance.
Even if the worst case could be computed accurately while maintaining the target frame rate, in the majority of cases the simulation would use but a fraction of the resources allocated to it.
This inefficiency comes at the cost of taking resources away from other sub-systems that could do with it, impacting the quality of the whole application.

---

There is, thankfully, a solution.
While Newton-Raphson is likely to have the best best-case performance - scenarios with nearly-circular orbits - it has appalling worst-case performance.
As soft real-time systems, games need to be optimised for the worst-case, such that framerates are consistently high.
We also want to make the most of our available resources on every frame, rather than leave most of them unused most of the time.
For this, we can turn to that old favourite of computer scientists: a binary search.

```c
Real M = meanAnomaly;
Real e = eccentricity;
Real E;
Real D;

Real F = sign(M);
M = abs(M) / (2.0 * pi);
M = (M - (int)M) * 2.0 * pi * F;
if (M < 0)
{
    M = M + 2.0 * pi;
}

F = 1.0;
if (M > pi)
{
    F = -1.0;
    M = 2.0 * pi - M;
}

E = pi * 0.5;
D = pi * 0.25;
for (int J = 0; J < iterations; J++)
{
    Real M1 = E - e * sin(E);
    Real val = M - M1;
    E = E + D * sign(val);
    D *= 0.5;
}

E *= F;

return E;
```

Eccentric anomaly determination from _Meeus, J. Astronomical Algorithms. Second Edition, Willmann-Bell, Inc. 1991._
Translated to pseudo-C from BASIC source.

---

This approach, while obtuse and archaic, has the property of running in a predictable, constant amount of time, and computes a result to a particular level of precision determined by the number of iterations.
On average, it performs worse than Newton-Raphson, taking more runtime to arrive at a result of similar accuracy.
However, in the worst-case, it performs just as well in terms of accuracy and runtime as the best-case, which makes it suitable for real-time applications in a way that Newton-Raphson is not.
Instead of having to pad the physics budget to account for a worst-case scenario, we can make much more efficient use of a smaller budget, despite the algorithm almost always costing more.

---

I don't know how Kerbal Space Program computes positions on elliptical orbits - if you do, [please get in touch!](mailto:me@ktyl.dev) - but I would be surprised if they used Newton-Raphson.
In practice, Kepler's laws are limited in their utility for complex, real-life astrodynamical systems, as they don't take into account tidal effects, irregular gravitational fields, atmospheric drag, solar pressure or general relativity.
They also cannot model Lagrange points, which are fundamental to many modern spacecraft such as JWST, or Weak Stability Boundary transfers as used by modern robotic lunar probes such as SLIM.
Nonetheless, game developers are a resourceful bunch, so I'm still holding my breath for a game about the [Interplanetary Transport Network](https://en.wikipedia.org/wiki/Interplanetary_Transport_Network).