25
$\begingroup$

This might come of as a millionth-time award winner beginner question, but since I'm not familiar with either correct terms or designations, I'll put this question in my own language, since I'm not really sure what to search for in previously answered questions.

I am building a game where I simulate planetary orbits. These orbits are static and will not change, as expected. Up until now I have calculated these orbits with step-by-step iteration with gravity vectors, but the scale have grown beyond what is possible in real-time.

My question is therefore how I calculate the planet's position at a given time $t$:

Known variables:

  • planet position at time 0
  • planet speed at time 0
  • planet direction at time 0

What I wish to know:

  • planet position at time $t$

I've probably seen some fancy answers to this questions out there, but they are not very intuitive for someone with limited knowledge in astrophysics like me. Thank you for your patience with me!

$\endgroup$
38
$\begingroup$

Statement of the Problem

The problem you want to solve is called the Kepler problem. In your formulation of the problem, you're starting out with the Cartesian orbital state vectors (also called Cartesian elements): that is, the initial position and velocity.

As you have discovered, the only way to propagate the Cartesian elements forward in time is by numerical integration. This works OK, but it can be slow if you want high accuracy, and there are some numerical problems (errors caused by rounding[slowly accumulate, and many integrators cause energy drift). You can get around some of these problems by using a higher-order integrator (Runge-Kutta is a popular one) which allows you to take larger steps for the same level of accuracy, or get better accuracy for the same step size. However, this is kind of overkill for simple simulation.

If your simulation can be treated as a two-body problem, then things simplify dramatically. The two-body problem is a good simplification if the simulation objects are mainly affected by a single, large object. For example, the Earth traveling around the Sun or a spacecraft traveling in a low Earth orbit are well-modeled as a two body problem; however, a spacecraft traveling from the Earth to the Moon is not (more on that later).

Since you're trying to model the positions of the planets with medium accuracy, reduction to the two-body problem should work for you.

Definition of Terms

The traditional solution to the two-body problem involves a different way of representing the position of the orbiting body, called the Keplerian orbital elements (also just called orbital elements). Instead of specifying position and speed, they specify six different parameters of the orbit (if you just want to get to the code, you can skip this part):

  • Semi-major axis, $a$: Half the length of the orbit, analogous to radius. The energy and period of the orbit depend only on $a$. The semi-latus rectum $\ell$, the "width" of the orbit, can be a better choice for orbits that are close to parabolic (like for asteroids) or that change from elliptical to hyperbolic (like interplanetary spacecraft). The two are related by $\ell=a(1-e^2)$.

  • Eccentricity, $e$: The "pointiness* of the orbit. Ranges from $e=0$ for a perfectly circular orbit, to $e=1$ for a parabolic orbit, to $e>1$ for hyperbolic orbits. Mercury is the most eccentric planet with $e \approx 0.2$. Earth-orbiting spacecraft usually have $e<0.01$.

    Aside: From $e$ and $a$ we can determine the highest and lowest points in the orbit, the apoapsis and periapsis (together apsides): $$ r_a = a(1+e) \\ r_p = a(1-e) $$ The naming of these points is a little funny: apoapsis and periapsis are the generic terms, but orbits around particular bodies have specific terms: a spacecraft around the Earth has an apogee and perigee, while the Earth (in orbit around the Sun) has an aphelion and perihelion.

The two parameters $a$ and $e$ are enough to determine the shape of the orbit. The next three parameters define the orientation of the orbit relative to a coordinate system consisting of a reference plane, and a reference direction (parallel to the plane).

For almost all orbits in the solar system, the coordinate system used is the ecliptic coordinate system. The reference plane is the ecliptic plane, the plane of the Earth's orbit around the Sun. The reference direction is the vernal equinox point, the direction from the Earth to the Sun at the moment of the spring equinox. Since both these references drift slowly over time, we must specify a particular time at which these references are defined, called the epoch. The most common is J2000, noon on January 1, 2000 (UTC).

Earth-centered orbits often use the equatorial coordinate system, whose reference plane is the equator of the Earth. The situation with the epoch is a little complicated, so I won't get into it here.

  • Inclination, $i$: the angle between the plane of the orbit and the reference plane. Inclination between 90 and 180 degrees refers to a retrograde orbit, one that orbits "backwards" from the usual direction.

  • Longitude of the ascending node, $\Omega$: the ascending node is where the orbit crosses from below the reference plane to above it. (It's at the intersection between the orbital plane and the reference plane) $\Omega$ is the angle between this point and the reference direction, measured counterclockwise.

  • Argument of periapsis, $\omega$: the angle between the ascending node and periapsis (the lowest point in the orbit). For orbits with very low inclination where the location of the ascending node is difficult to determine (since it's the intersection between two almost parallel planes), we instead use the longitude of periapsis $\varpi = \Omega + \omega$.

The sixth parameter defines the position of the object in its orbit. There are a couple different choices, but the most common is:

  • Mean anomaly, $M$: an "imaginary" angle that is zero at periapsis and increases at a constant rate of 360 degrees per orbit.

The rate at which $M$ changes is called the mean motion, $n$, equal to $2\pi/T$. Usually you have a measurement of $M$ at a particular epoch $t_0$, called (unsurprisingly) the mean anomaly at epoch, $M_0$.

Just like the argument of periapsis, for low-inclination orbits we use a related value, the mean longitude, $L=\varpi + M$.

The actual angle between the orbiting body and periapsis is called the true anomaly, $\nu$. This is the angle we need to compute the position of the body. Unfortunately there is no way to directly compute $\nu$ from $M$. Instead we first solve for the eccentric anomaly $E$:

$$ M = E - e \sin E $$

This is called Kepler's equation, and it cannot be solved analytically. Once we have $E$ though, there is a relatively simple expression for $\nu$.

Computing Position from Orbital Elements

We'll perform this computation in three steps: first, we'll solve Kepler's equation. Second, we'll compute the 2d position of the body in the orbital plane. Lastly, we'll rotate our 2d position into 3d coordinates. I'll give some "pseudocode" in Javascript for most of these tasks.

I'll assume that you're using a set of elements like these from JPL's web site. These use $L$ and $\varpi$ instead of $M$ and $\omega$. The table gives two values for each of the elements; the second is the time derivative. If you use the values in this table you should use the derivatives as well.

Calculate the time $t$ in centuries from J2000:

// month is zero-indexed, so 0 is January
var tMillisFromJ2000 = Date.now() - Date.UTC(2000, 0, 1, 12, 0, 0);
var tCenturiesFromJ2000 = tMillisFromJ2000 / (1000*60*60*24*365.25*100);

Now we calculate the current values of each of the orbital parameters. For example, the semimajor axis of Earth, using the values from Table 1 (valid from 1800–2500):

// a0 = 1.00000261; adot = 0.00000562
var a = a0 + adot * tCenturiesFromJ2000;

(Note that the values are actually given for "EM Barycenter," the center-of-mass of the Earth-Moon system. The Earth is around 4600 kilometers from the barycenter in the opposite direction from the Moon. If you want to correct this inaccuracy you'll need to simulate the motion of the Moon as well, but that's probably overkill.)

Table 2a gives elements that are accurate from 3000 BC to 3000 AD; however, if you use the elements from table 2a, you must supplement them with corrections to $L$ from Table 2b! For example, here is computing the longitude of Saturn:

// L0 = 34.33479152; Ldot = 3034.90371757
// b = -0.00012452          
// c =  0.06064060
// s = -0.35635438
// f = 38.35125000
var L = L0 + Ldot * tCenturiesFromJ2000
           + b * Math.pow(tCenturiesFromJ2000, 2)
           + c * Math.cos(f * tCenturiesFromJ2000)
           + s * Math.sin(f * tCenturiesFromJ2000);

We don't need to explicitly compute the mean motion and add it to $L$, since both tables include it in $\dot L$.

Now we're ready to compute $M$ and $\omega$ (w):

var M = L - p \\ p is the longitude of periapsis
var w = p - W \\ W is the longitude of the ascending node

On to step 2: we need to solve the Kepler equation:

$$ M = E - e \sin E $$

We can solve this numerically using Newton's method. Solving the Kepler equation is equivalent to finding the roots of $f(E) = E - e \sin E - M$. Given $E_i$, an estimate of $E$, we can use Newton's method to find a better estimate:

$$ E_{i+1} = E_i - f(E_i) / f'(E_i) \\ f'(E) = 1 - e \cos E $$

Since the nonlinear part $e \sin E$ is very small, we can start with the estimate $E=M$. Our code looks something like this:

E = M;
while(true) {
  var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E));
  E -= dE;
  if( Math.abs(dE) < 1e-6 ) break;
}

Now there are two ways to compute the position from the eccentric anomaly. We can first compute the true anomaly and radius (the position of the object in polar coordinates), and then convert to rectangular coordinates; however, if we apply a bit of geometry we can instead compute the coordinates directly from $E$:

var P = a * (Math.cos(E) - e);
var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2));

(P and Q form a 2d coordinate system in the plane of the orbit, with +P pointing towards periapsis.)

Finally, we can rotate these coordinates into the full 3d coordinate system:

// rotate by argument of periapsis
var x = Math.cos(w) * P - Math.sin(w) * Q;
var y = Math.sin(w) * P + Math.cos(w) * Q;
// rotate by inclination
var z = Math.sin(i) * x;
    x = Math.cos(i) * x;
// rotate by longitude of ascending node
var xtemp = x;
x = Math.cos(W) * xtemp - Math.sin(W) * y;
y = Math.sin(W) * xtemp + Math.cos(W) * y;

(x, y, and z will be in units of AU.)

And you're done!


A few tips:

  • If you want to calculate the velocity as well, you can do it at the same time as you calculate $P$ and $Q$, then rotate it in the same way. $$ \dot M = n = \dot L \\ \dot M = \dot E - e (\cos E) \dot E \\ \dot E = \dot M / (1 - e \cos E) \\ \dot P = -a (\sin E) \dot E \qquad \dot Q = a (\cos E) \dot E \sqrt{1 - e^2} $$ Note I don't include any of the derivatives (except $\dot L$) in this calculation, since they don't affect the outcome much. You could code this as:

    var vP = - a * Math.sin(E) * Ldot / (1 - e * Math.cos(E));
    var vQ = a * Math.cos(E) * Math.sqrt(1 - e*e) * Ldot / (1 - e * Math.cos(E));
    

    Note that the velocities will be in AU per century.

  • If you are updating the positions very frequently, you could use the previous value of $E$ to seed Newton's method, and do a fixed number of iterations (probably just one would suffice). Note however that you need to keep that value of $E$ local to each object!

  • You can also just use a fixed number of iterations for the initial solution. Even for $e=0.2$, after three iterations the error in $E$ is only about $10^{-13}$, and after four iterations the error is smaller than the rounding error of an IEEE double up to $e=0.42$.

If you want more information you can search online, but if you're really interested you should read an introductory text on orbital mechanics. I personally recommend Fundamentals of Astrodynamics by Bate, Mueller, and White (pdf). My dad used this book back when he was in college, and I found it to be more readable than my college textbook. You'd be interested in Chapter 4, Position and Velocity as a Function of Time.

$\endgroup$
  • 1
    $\begingroup$ Here are the orbital elements of the planets to get you started. $\endgroup$ – 2012rcampion Apr 26 '15 at 15:23
  • $\begingroup$ It should probably be mentioned that the system above is for "two body motion": the planets don't affect each other. if the original poster's program used to be including these effects, then the results will be somewhat different due to this issue. $\endgroup$ – Matt Jessick Feb 4 '16 at 0:25
  • $\begingroup$ @Matt good point, I'll add that in. $\endgroup$ – 2012rcampion Feb 4 '16 at 1:20
  • $\begingroup$ Please, can you review my/your code? I translated it to C# and it didn't work. Info about the problem here: space.stackexchange.com/questions/21458/… Full, updated code here: pastebin.com/ydse66XY $\endgroup$ – vistaero May 13 '17 at 17:38
  • 1
    $\begingroup$ @vistaero You need to convert the angles to radians before using Keper's equation; otherwise you need to add a factor of $180/\pi$ to the $e \sin E$ term and its derivative. $\endgroup$ – 2012rcampion May 14 '17 at 0:42
8
$\begingroup$

Since its just a game, would you be happy with circular orbits and the planets' orbits only affected by the central body? In that case, the propagation is quite simple. In the plane of the orbit with the central body at (0,0), the position as a function of time is:

$$x(t)=a\cos\left({2\pi \left(t-t_0\right)\over T}\right)$$

$$y(t)=a\sin\left({2\pi \left(t-t_0\right)\over T}\right)$$

where $a$ is the semi-major axis, or really just the orbit radius in this case, $T$ is the orbit period, and $t_0$ determines the phasing of the orbit, where at $t=t_0$, the planet is on the x-axis on the positive side.

To make the orbits of the different planets consistent with each other, you just need to define the $GM$ of the central body, which we will call $\mu$. Then for any orbit radius $a$, the orbit period is related to $a$ by:

$$T=2\pi\sqrt{a^3\over\mu}$$

$\endgroup$
2
+150
$\begingroup$

While there has already been a high quality accepted answer for years, here is some additional background, some particularly helpful resources, and additional tips for first-time orbit propagation.


If you're not doing N-body physics, so the planets do not interact then you can use analytic solutions to the Kepler problem. Eventually you'll realize that you need to solve hyperbolic orbits at some point as well. That will lead you to universal variables formulations of solving the Kepler problem.

The best solutions to that are probably going to be Goodyear's method:

W. Goodyear, “Completely General Closed Form Solution for Coordinates and Partial Derivatives of the Two-Body Problem”, The Astronomical Journal, Vol. 70, No. 3, 1965, pp. 189–192 (or the NASA NTRS TD document on the same material)

Shepperd's method:

Shepperd, S.W. Celestial Mechanics (1985) 35: 129. https://doi.org/10.1007/BF01227666

Or Danby-Stumpff:

Danby, J.M.A. Celestial Mechanics (1987) 40: 303. https://doi.org/10.1007/BF01235847

There is some MATLAB code here which might be useful (and vastly more accessible), although random code snippets on matlabcentral are far from guaranteed to be bug free and it looks like this code may lack useful normalization of its inputs (generally you're going to want to normalize to the scale of your problem so that you do math in units where r0-bar = 1.0 and mu-bar = 1.0 and where v-bar = 1 is the velocity in a circular orbit at r0 or something like that).

If you are going to do N-body integration of planetary motion then I think you're going to have to use numerical integration. Runge-Kutta will violate conservation of Energy so you will likely want to use Symplectic Integration. The 4th order symplectic integrator in that article is not that difficult to code -- although that leaves you with the difficulty of guessing the correct timestep (again, normalization helps because a circular planetary orbit and circular LEO are the same problem just with different distance scales) and with interpolation of the interior points (and you need to watch out for Runge's phenomenon, but I haven't wrestled with that, so don't know which approach to take there).

If you're going to use Runge-Kutta then Dormand-Prince with dynamic step side and its 3rd order interpolant will be very convenient, and is what Matlab uses in its ode45 solver.

I would probably advise starting with the simplest runge-kutta implementation based on ease of coding, but if you're doing runge-kutta on every physics tick to advance it forward one step then that is pretty brutal and the errors will eventually add up, but you could prototype it that way. At some point you'll want to go to a system where you solve the problem for many time steps into the future, and then you use an interpolating function to pick off the solution at intermediate timesteps (which is the point of my mentioning Dormand-Prince and its interpolating function).

$\endgroup$
  • 1
    $\begingroup$ This is a great, and very helpful answer to an old question. Thanks for the excellent sources; I'm going to look these up in the next few days. $\endgroup$ – uhoh Oct 18 '18 at 8:24
  • 1
    $\begingroup$ Yeah, I've worked with KSP's orbit class for years, and know something about Principia is built. I would strongly suggest for anyone setting out to build a game, particularly if you don't just want to clone KSP but try to build something better, that these issues get thought about up front in a bit more detail. Looks like I was in a hurry yesterday and missed "these [planetary] orbits are static" so maybe all my yapping about N-body and Runge-Kutta methods was superfluous, but other people who find this question might care. And Goodyear's solution still beats KSP's Orbit class for speed. $\endgroup$ – lamont Oct 18 '18 at 17:34
  • 1
    $\begingroup$ (the other solutions after Goodyear's solution may also be even faster, I just don't have experience with those) $\endgroup$ – lamont Oct 21 '18 at 16:39
  • $\begingroup$ I'd love to read those, but they are paywalled. I'll check them out the next time I get to the library. $\endgroup$ – uhoh Oct 21 '18 at 16:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.