In the previous post, we saw how to do the math to get velocities and positions along an elliptical orbit. To translate this into game code though, there is still one problem we have to solve. So far we have spoken about apoapsis and periapsis as determining the shape of the ellipse that describes the orbit, but it doesn't uniquely identify the orbit itself. As you can see in the image below, the rendered orbits have the same periapsis and apoapsis, but they still do not coincide. The orbit has been rotated around the primary, as it were.
To uniquely identify orbits, we need one more parameter, namely where the orbit reaches periapsis. We talk about the longitude of periapsis $\varpi$ with respect to some reference direction. For our code, we will assume the positive x-axis as our reference direction.
We assume that the satellite is at periapsis at some instant $t_{0}$, and we now know that the relative angle with the x-axis at periapsis is $\varpi$. This means we have our starting point, and from there we can use Kepler's equations to calculate the rest.
Current position in elliptical orbit
In the previous post, we showed how to get the current true anomaly based on the elapsed time, the eccentricity, and mean anomaly of the orbit. As a reminder, the mean anomaly is the average change in the direction of the satellite around the primary, so it is just $2\pi$ (a full rotation) divided by the period. To find the period of an elliptical orbit, we look at our good old friend Kepler again. In particular, Kepler's Third Law says:
$$T = 2\pi \sqrt{\frac{a^3}{GM}}.$$
$T$ is the period in seconds, $a$ is the semi-major axis of our orbit as before, and $GM$ the gravitational constant multiplied with the mass of the primary. As with circular orbits before, as a game developer you have the option to choose these however you want to make things look as you want to. The mean anomaly is thus:
$$\mu = \frac{2\pi}{T} = \frac{2\pi}{2\pi \sqrt{\frac{a^3}{GM}}} = \frac{1}{\sqrt{\frac{a^3}{GM}}}$$
We now have a plan: whenever we create a new orbit in our game, with a primary, periapsis, apoapsis, $t_0$, and longitude of periapsis, we first need to calculate a few derivative properties: the eccentricity, the semi-major and semi-minor axes, the centre of the ellipse, and the mean anomaly. These values don't change. From that point onward, we use Kepler's equation to calculate the position ever frame.
The full code can be found below.
static class Constants {
public const double G = 1;
}
sealed class Body {
public Vector2 Position { get; }
public double Mass { get; }
}
sealed class Orbit {
private readonly Body primary;
private readonly double tPer;
private readonly double eccentricity;
private readonly double semiMajorAxis;
private readonly double semiMinorAxis;
private readonly Vector2 ellipseCenterOffset;
private readonly double meanAnomaly;
private readonly Matrix2 longitudeTransform;
public Orbit(
Body primary, double periapsis, double apoapsis, double tPer, double longitudePer) {
this.primary = primary;
this.tPer = tPer;
eccentricity = (apoapsis - periapsis) / (apoapsis + periapsis);
semiMajorAxis = 0.5 * (periapsis + apoapsis);
semiMinorAxis = semiMajorAxis * Math.Sqrt(1 - eccentricity * eccentricity);
ellipseCenterOffset = (periapsis - semiMajorAxis) * Vector2.UnitX;
meanAnomaly = 1 / Math.Sqrt(
semiMajorAxis * semiMajorAxis * semiMajorAxis / (Constants.G * primary.Mass));
longitudeTransform = Matrix2.CreateRotation(longitudePer);
}
public Vector2 PositionAtTime(double currentTime)
{
// Calculate the true anomaly using Kepler's equation.
var trueAnomaly = trueAnomalyAtTime(currentTime);
// We calculate the distance vector from the primary to the satellite
// under the assumption that longitude of periapsis is 0.
var preRotatedOffset = ellipseCenterOffset + new Vector2(
semiMajorAxis * Math.Cos(trueAnomaly),
semiMinorAxis * Math.Sin(trueAnomaly));
// We then rotate the point around based on the longitude of periapsis.
var rotatedOffset = longitudeTransform.Transform(preRotatedOffset);
// Finally, we translate this offset to be in the space of the primary.
return primary.Position + rotatedOffset;
}
private double trueAnomalyAtTime(double currentTime)
{
const numIterations = 30; // Change depending on expected eccentricities.
// Pre-calculate the current mean motion because it doesn't change.
var M = (currentTime - tPer) * meanAnomaly;
// Make an initial guess based on the eccentricity of the orbit.
var E = eccentricity > 0.8 ? Math.PI : M;
for (int i = 0; i < numIterations; i++)
{
var numerator = eccentricity * Math.Sin(E) + M - E;
var denominator = eccentricity * Math.Cos(E) - 1;
E = E - numerator / denominator;
}
return E;
}
}
That is a lot of code, but it is all derived directly from the equations we've seen over the past few posts. There is one new trick I have applied to deal with the longitude of periapsis though. Calculating offsets is hard when your ellipse is not axis-aligned. That is why for all the calculations, I just pretend our primary is the origin, and the periapsis lies along the positive x-axis. When we calculated the difference vector from the primary to the satellite in this normalized space, we transform it by first rotating around the origin, and then calculating the final position by adding the rotated difference vector to the primary's position.
The third dimension
We spent five posts getting to the point of being able to describe an elliptical orbit in code. Surely, doing this in three dimensions will be much harder! Nothing less would be the case. All the equations still hold up in three dimensions. As a matter of fact, we have assumed that the physical rules of a 3D world hold up, and that we are merely working in a two-dimensional projection of the 3D world. This is also exactly why the step from 2D to 3D is relatively straightforward: we just have to put our code to work in a specific plane in that 3D world.
In orbital mechanics, we don't only have a reference for the longitude of periapsis, we also have a reference plane. In games we could easily pick the xy-plane, in the real world it is often the plane that is set out by the primary's equator. Our orbit will either coincide with the plane, or it will pass the plane through two points. The angle the orbit makes with the reference plane is called the inclination. The points in which the satellite passes through the reference plane are called the ascending node (when the orbit goes from below the plane to above it) and descending node (when the orbit goes in the other direction).
Earlier we used the longitude of periapsis to determine the orientation of the orbit, but in practice, we don't actually use that parameter. Instead, we use something called the longitude of the ascending node $\Omega$. Instead of measuring the angle between the reference direction and the periapsis, we measure the angle between the reference direction and the ascending node. To then fix where the periapsis is, we express the position of the periapsis as the argument of periapsis $\omega$ (yes they're both omegas, I didn't choose the symbols). The argument of periapsis is the angle between the ascending node and the periapsis measured in the orbital plane. The reason it is chosen this way is that with third dimensions, the periapsis no longer necessarily lies in the reference plane, meaning that calculating the angle between a reference direction and the periapsis becomes an undefined problem.
To implement the longitude of periapsis in our code, we used a transformation that kept the math simple by moving the primary to the origin, and the periapsis on the positive x-axis. We can do the same thing for inclination. Remember that transformations are always applied "backwards" or right to left, so the steps in the code become:
- Rotate
rotatedOffsetaround the z-axis (so it rotates in the reference plane xy) by the argument of periapsis. Our orbit now represents how the orbit looks with the orbital plane as reference plane, and the ascending node at the positive x-axis. - Rotate
rotatedOffsetaround the x-axis by the inclination. We now have our actual orbit, but the ascending node is not yet in the right position. - Rotate
rotatedOffsetaround the z-axis, since the xy-plane now the actual reference plane, by the longitude of the ascending node.
Since these three rotations never change, we can precalculate the resulting matrix, and just translate our result once. That makes our final code:
static class Constants {
public const double G = 1;
}
sealed class Body {
public Vector2 Position { get; }
public double Mass { get; }
}
sealed class Orbit {
private readonly Body primary;
private readonly double tPer;
private readonly double eccentricity;
private readonly double semiMajorAxis;
private readonly double semiMinorAxis;
private readonly Vector2 ellipseCenterOffset;
private readonly double meanAnomaly;
private readonly Matrix2 orbitTransform;
public Orbit(
Body primary, double periapsis, double apoapsis,
double tPer, double inclination, double longitudeAscNode, double argPeriapsis) {
this.primary = primary;
this.tPer = tPer;
eccentricity = (apoapsis - periapsis) / (apoapsis + periapsis);
semiMajorAxis = 0.5 * (periapsis + apoapsis);
semiMinorAxis = semiMajorAxis * Math.Sqrt(1 - eccentricity * eccentricity);
ellipseCenterOffset = (periapsis - semiMajorAxis) * Vector2.UnitX;
meanAnomaly = 1 / Math.Sqrt(
semiMajorAxis * semiMajorAxis * semiMajorAxis / (Constants.G * primary.Mass));
orbitTransform = createOrbitalTransformMatrix(
inclination, longitudeAscNode, argPeriapsis);
}
private Matrix3 createOrbitalTransformMatrix(
double inclination, double longitudeAscNode, double argPeriapsis) {
// Rotate the periapsis in position assuming the ascending node is the
// positive x-axis.
var periapsisTransform = Matrix3.CreateAxisRotation(Vector3.UnitZ, argPeriapsis);
// Rotate the orbit around the x-axis to give it its inclination.
// We use the negative rotation here to make sure the x-axis is the
// ascending node.
var inclinationTransform = Matrix3.CreateAxisRotation(Vector3.UnitX, -inclination)
// Finally, rotate the orbit so that the ascending node is in the right
// position as well.
var ascNodeTransform = Matrix3.CreateAxisRotation(Vector3.UnitZ, longitudeAscNode);
// Transformations are applied right to left, so make sure we multiply
// them in that order.
return ascNodeTransform * inclinationTransform * periapsisTransform;
}
public Vector2 PositionAtTime(double currentTime)
{
// Calculate the true anomaly using Kepler's equation.
var trueAnomaly = trueAnomalyAtTime(currentTime);
// We calculate the distance vector from the primary to the satellite
// under the assumption that longitude of periapsis is 0.
// Note that the z-coordinate is 0, since the orbit is in the xy-plane.
var preRotatedOffset = ellipseCenterOffset + new Vector3(
semiMajorAxis * Math.Cos(trueAnomaly),
semiMinorAxis * Math.Sin(trueAnomaly),
0);
// We then rotate the point in 3D space to apply inclination and the
// correct orientation of the periapsis and ascending node.
var rotatedOffset = orbitTransform.Transform(preRotatedOffset);
// Finally, we translate this offset to be in the space of the primary.
return primary.Position + rotatedOffset;
}
private double trueAnomalyAtTime(double currentTime)
{
const numIterations = 30; // Change depending on expected eccentricities.
// Pre-calculate the current mean motion because it doesn't change.
var M = (currentTime - tPer) * meanAnomaly;
// Make an initial guess based on the eccentricity of the orbit.
var E = eccentricity > 0.8 ? Math.PI : M;
for (int i = 0; i < numIterations; i++)
{
var numerator = eccentricity * Math.Sin(E) + M - E;
var denominator = eccentricity * Math.Cos(E) - 1;
E = E - numerator / denominator;
}
return E;
}
}
For bonus points, since all the transformations are rotations, you could switch to using a quaternion instead to speed up that step, but that part is left as an exercise to the reader.
Retrograde orbits
For the last few posts, we've been silently assuming that the orbits move counter-clockwise. If our primary also rotates around its axis counter-clockwise, we call these prograde orbits (i.e. the direction of the rotation of the primary matches the direction of the orbit). Prograde orbits are much more common than orbits in the opposite direction: retrograde orbits (i.e. orbits that move in the opposite direction as that of rotation of the primary). However, now that we have inclinations, retrograde orbits are a breeze to implement. To make an orbit traverse backwards, we just rotate the entire orbit by 180 degrees around some point.
Conclusion
It took us several blog posts to get to this point, but we are now able to implement orbits using a somewhat analytical approach. We started out by seeing how integration methods can cause our orbital bodies to drift, and eventually just fly off the screen due to accumulated error. While recalculating the position of a satellite every frame is more work, starting from scratch each frame means we don't carry over our previous errors. The satellite may still not be in exactly the right position, but at least the error will not increase in magnitude as our simulation runs longer.
This solution is great, and maybe even necessary, for games that rely on stable orbits, but the solution doesn't scale. Moving from a two-body orbital system to a three- or even n-body orbital system is nearly impossible, and it is likely you will need to fall back on the physical simulation by frame method in those cases.
In case you want to read any of the past posts on this topic, below is a list of all posts in this series.