Kepler Orbits And Exoplanets

Lately, I've been working on finding exoplanet orbits based on their light data, and so I ended up with some code that I decided to turn into a little gadget that renders Keplerian orbits using Rust and WebAssembly.

Why did I decide to do it when NASA already has a much better one on their website? I'll tell you:

This post covers three things:

1. How do Keplerian orbits work?
2. How can you use this knowledge to find Exoplanets?
3. How don't Keplerian orbits work: wobble, drift, etc.

I. Keplerian Orbits

Below you can see what is meant by Keplerian orbits. It's a tiny model of the Solar System with the planets orbiting the Sun. You can drag the scene below and tweak the orbital parameters with the widget below. You can also zoom in and out with your mouse wheel.

I've written solutions to Keplerian orbits three times in my life, and the first two times I got it wrong. It turns out, turning 16th century equations into javascript that spins the planets the right way around is easy to mess up.

It doesn't help that many of the examples you find online are straight up wrong, and it's hard to tell which ones. The problem is that when you look up the solution to Kepler Orbits, you typically find something like this concise document.

Not only is this needlessly complicated, it's also easy to get wrong for programmers, because computers and math don't mix as well as you might think, and floats in C++ are, like, whatever.

Here's some ways to mess up:

• Things move in the wrong direction, while other things are moving in the right direction.
• Things move fast when they should be moving slowly, and vice-versa.
• Precision errors will creep up around $\tan$ values, which tend to infinity with every period.
• One of the 50 $\sin$ and $\cos$ calls will have the wrong sign or use $\varpi$ instead of $\omega$, resulting in complete and inexplicable chaos.

Anyway...

It need not be this way. You can get a solution that's easy to understand and obviously correct. We just break it up into three steps:

1. Draw an ellipse
2. Find the planet based on its starting position
3. Rotate the orbit in 3D using two angles

For most of this, you only need high school math. (Solving Kepler's Equation technically requires knowledge of Newton's Method, but just guessing is good enough. Rotation into 3D is made a whole lot easier if you know quaternions, but it's doable with just a lot of trig, too.)

The rest of this article will be easy to read if you've ever made a video game in Unity or Unreal, or if you were good at trig in 9th grade.

Orbital Elements

To describe the orbit of a planet around the Sun (or a satellite around the Earth), we require at least six parameters. There are various combinations to choose from, but the two we'll look at today are Keplerian elements.

Keplerian elements, at first glance, might look esoteric:

1. $a$ Semi-major axis - half of the width of the ellipse.
2. $e$ or $\varepsilon$ Eccentricity – half-open interval $[0, 1)$ measuring deviation from the perfect circle, 0 being a perfect circle and 1 being a parabola.
3. $M$ Mean anomaly – a longitudinal angle that measures where the planet would be when seen from the Sun, if the orbit were perfectly circular.
   • Alternatively, some sources (like NASA) give an angle called $L$ mean longitude, which is just $M + \varpi$.
4. $\varpi$ Longitude of the perihelion – the angle that measures where the planet makes its closest approach to the Sun, called perihelion¹.
5. $\Omega$ Longitude of the ascending node – the angle that measures where the planet's orbital plane intersects with the orbital plane of the Earth, called the ecliptic².
6. $I$ Inclination – the angle between the orbital plane and the ecliptic.

This is a lot of parameters, but, happily, we can deal with them in turns. First, we can deal with the first two, which just specify an ellipse.

Step 1: Just draw an ellipse?

To draw an orbit in 2 dimensions, you just draw an ellipse. For this you need to know only the eccentricity $e$ and the semi-major axis $a$. You also need to know the high-school equation for ellipse:

\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \quad \text{where} \quad b = a\sqrt{1 - e^2}

Here's what one looks like:

Figure 4: An elliptical orbit of two parameters.

We need to rotate it using one of the angular parameters ($\varpi$ - Longitude of the perihelion).

Figure 5: An elliptical orbit of three parameters.

Step 2: Where is the planet?

Now we need to find the planet. For this, we need to know only an angle of where it is in its orbit, such as the Mean anomaly $M$. Unfortunately, $M$ would only give us the real angle (called the true anomaly $\nu$) if the orbit were a circle. Since it's not a circle, we have to do some algebra:

The mean anomaly angle $M$ is related to an intermediate result $E$ by this equation.

M = E - e \sin(E)

Once we know $E$, we can compute the true angle $\nu$ simply:

\tan\frac{\nu}{2} = \sqrt{\frac{1+e}{1-e}} \tan\frac{E}{2}

The problem, then, is figuring the value of $E$. The thing is, there is no algebraic, closed form way to rearrange the equation and solve for $E$. We, therefore, do the next best thing and just guess.

We start by guessing that $E = M$. It doesn't, but having an estimate $E_0$ lets us test it, find how much we're off by, and then guess again. Formally:

$$ E_{n+1} = E_n - \frac{E_n - e \sin(E_n) - M}{1 - e \cos(E_n)} $$

If you're into that sort of thing, here's a bit of code I wrote in Rust that automates this guessing process:

// Convert mean anomaly (M) to eccentric anomaly (E). Both angles in radians.
//
// This implements a numeric solution using Newton's method.
pub fn eccentric_anomaly<T>(mean_anomaly: T, eccentricity: T) -> T
where
    T: num_traits::float::Float + From<f32> + num_traits::float::FloatConst,
{
    // M = E - e × sin(E)
    // - M: mean anomaly
    // - E: eccentric anomaly
    // - e: eccentricity
    //
    // Using Newton's method, minimize the error function:
    //
    // f(E) = M - (E - e × sin(E))
    //
    // First-order derivative:
    //
    // f'(E) = e × cos(E) - 1
    let tolerance = Into::<T>::into(1e-6);
    let epsilon = Into::<T>::into(1e-14);
    let mut eccentric_anomaly = mean_anomaly;
    for _ in 0..100 {
        let dy = eccentricity * T::cos(eccentric_anomaly) - T::one();
        if T::abs(dy) < epsilon {
            return eccentric_anomaly;
        }
        let y = mean_anomaly - (eccentric_anomaly - eccentricity * T::sin(eccentric_anomaly));
        let est = eccentric_anomaly - y / dy;
        if T::abs(est - eccentric_anomaly) <= tolerance {
            return eccentric_anomaly;
        }
        eccentric_anomaly = est;
    }
    eccentric_anomaly
}

The code guesses up to a hundred times, but when people did this on paper, 2-3 guesses were generally good enough to find the planet in the sky with a telescope.

Figure 6: An elliptical orbit showing the position of a planet.

Step 3: Make it 3D

Now we need to incline the orbit out of the plane it's in and into 3D space. For this, we need two angles: one that tells us how much to incline it (0 to 90 degrees - the Inclination $I$) and another one that tells us the axis: $\Omega$ Longitude of the ascending node.

This step involves a lot of trig, but if you're using a 3D library and have existing implementations of Quaternions, it's surprisingly simple to code:

pub fn orbital_to_ecliptic<T>(elements: Elements<T>) -> nalgebra::UnitQuaternion<T>
where
    T: num_traits::float::Float + From<f32> + nalgebra::RealField,
{
    let polar_axis = Vector3::z_axis();
    let periapsis_rot = nalgebra::UnitQuaternion::from_axis_angle(&polar_axis, elements.periapsis);

    let inclination_axis =
        nalgebra::UnitQuaternion::from_axis_angle(&polar_axis, elements.ascending_node)
            * Vector3::x_axis();
    periapsis_rot
        * nalgebra::UnitQuaternion::from_axis_angle(&inclination_axis, elements.inclination)
}

Here's a single orbit with an inclination of about 45 degrees. Like before, you can drag to rotate the view and see how the orbit is tilted relative to the ecliptic plane. You can also tweak the parameters to get a better idea of what each one does.

Getting the orbital parameters of each planet.

NASA is a good source, as with many things.

Using this knowledge to find exoplanets

Read in Part II.

Exceptions to Kepler's laws

Read in Part III.

References

Figure 2: Lucas Snyder: "Diagram illustrating and explaining various terms in relation to Orbits of Celestial bodies." Accessible on Wikipedia. Licensed under Creative Commons Attribution-Share Alike 3.0 Unported.