Determining orbital position at a future point in time - Space Exploration Stack Exchange most recent 30 from space.stackexchange.com 2019-08-24T22:32:39Z https://space.stackexchange.com/feeds/question/8911 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://space.stackexchange.com/q/8911 26 Determining orbital position at a future point in time Henrik https://space.stackexchange.com/users/9021 2015-04-26T07:54:45Z 2019-08-15T08:41:43Z <p>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.</p> <p>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. </p> <h2>My question is therefore how I calculate the planet's position at a given time $t$:</h2> <p>Known variables:</p> <ul> <li>planet position at time 0</li> <li>planet speed at time 0</li> <li>planet direction at time 0</li> </ul> <hr> <p>What I wish to know:</p> <ul> <li>planet position at time $t$</li> </ul> <hr> <p>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!</p> https://space.stackexchange.com/questions/8911/-/8915#8915 41 Answer by 2012rcampion for Determining orbital position at a future point in time 2012rcampion https://space.stackexchange.com/users/8996 2015-04-26T15:20:25Z 2016-02-04T05:06:33Z <h2>Statement of the Problem</h2> <p>The problem you want to solve is called the <strong>Kepler problem</strong>. In your formulation of the problem, you're starting out with the <a href="http://en.wikipedia.org/wiki/Orbital_state_vectors"><strong>Cartesian orbital state vectors</strong></a> (also called <em>Cartesian elements</em>): that is, the initial position and velocity.</p> <p>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 <a href="https://en.wikipedia.org/wiki/Energy_drift">energy drift</a>). You can get around some of these problems by using a higher-order integrator (<a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods">Runge-Kutta</a> 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.</p> <p>If your simulation can be treated as a <a href="https://en.wikipedia.org/wiki/Two-body_problem"><strong>two-body problem</strong></a>, 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).</p> <p>Since you're trying to model the positions of the planets with medium accuracy, reduction to the two-body problem should work for you.</p> <h2>Definition of Terms</h2> <p>The traditional solution to the two-body problem involves a different way of representing the position of the orbiting body, called the <a href="http://en.wikipedia.org/wiki/Orbital_elements#Keplerian_elements"><strong>Keplerian orbital elements</strong></a> (also just called <em>orbital elements</em>). 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):</p> <ul> <li><p><a href="https://en.wikipedia.org/wiki/Semi-major_axis"><strong>Semi-major axis, $a$</strong></a>: Half the length of the orbit, analogous to radius. The energy and period of the orbit depend only on $a$. The <a href="https://en.wikipedia.org/wiki/Conic_section#Conic_parameters"><em>semi-latus rectum</em> $\ell$</a>, 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)$.</p></li> <li><p><a href="https://en.wikipedia.org/wiki/Orbital_eccentricity"><strong>Eccentricity, $e$</strong></a>: The "pointiness* of the orbit. Ranges from $e=0$ for a perfectly circular orbit, to $e=1$ for a parabolic orbit, to $e&gt;1$ for hyperbolic orbits. Mercury is the most eccentric planet with $e \approx 0.2$. Earth-orbiting spacecraft usually have $e&lt;0.01$.</p> <p>Aside: From $e$ and $a$ we can determine the highest and lowest points in the orbit, the <em>apoapsis</em> and <em>periapsis</em> (together <a href="https://en.wikipedia.org/wiki/Apsis"><em>apsides</em></a>): $$r_a = a(1+e) \\ r_p = a(1-e)$$ <sub>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 <em>apogee</em> and <em>perigee</em>, while the Earth (in orbit around the Sun) has an <em>aphelion</em> and <em>perihelion</em>.</sub></p></li> </ul> <p>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).</p> <p>For almost all orbits in the solar system, the coordinate system used is the <a href="https://en.wikipedia.org/wiki/Ecliptic_coordinate_system"><em>ecliptic coordinate system</em></a>. The reference plane is the <a href="https://en.wikipedia.org/wiki/Ecliptic"><em>ecliptic plane</em></a>, the plane of the Earth's orbit around the Sun. The reference direction is the <a href="https://en.wikipedia.org/wiki/Equinox#Celestial_coordinate_systems"><em>vernal equinox point</em></a>, 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 <a href="https://en.wikipedia.org/wiki/Epoch_(astronomy)"><strong>epoch</strong></a>. The most common is <a href="https://en.wikipedia.org/wiki/Epoch_(astronomy)#Julian_years_and_J2000">J2000</a>, noon on January 1, 2000 (UTC).</p> <p>Earth-centered orbits often use the <a href="https://en.wikipedia.org/wiki/Equatorial_coordinate_system"><em>equatorial coordinate system</em></a>, 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.</p> <ul> <li><p><a href="https://en.wikipedia.org/wiki/Orbital_inclination"><strong>Inclination, $i$</strong></a>: the angle between the plane of the orbit and the reference plane. Inclination between 90 and 180 degrees refers to a <em>retrograde</em> orbit, one that orbits "backwards" from the usual direction.</p></li> <li><p><a href="https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node"><strong>Longitude of the ascending node, $\Omega$</strong></a>: 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.</p></li> <li><p><a href="https://en.wikipedia.org/wiki/Argument_of_periapsis"><strong>Argument of periapsis, $\omega$</strong></a>: 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 <em>longitude of periapsis</em> $\varpi = \Omega + \omega$.</p></li> </ul> <p>The sixth parameter defines the position of the object in its orbit. There are a couple different choices, but the most common is:</p> <ul> <li><a href="https://en.wikipedia.org/wiki/Mean_anomaly"><strong>Mean anomaly, $M$</strong></a>: an "imaginary" angle that is zero at periapsis and increases at a constant rate of 360 degrees per orbit.</li> </ul> <p>The rate at which $M$ changes is called the <a href="https://en.wikipedia.org/wiki/Mean_motion"><em>mean motion</em>, $n$</a>, equal to $2\pi/T$. Usually you have a measurement of $M$ at a particular epoch $t_0$, called (unsurprisingly) the <em>mean anomaly at epoch</em>, $M_0$.</p> <p>Just like the argument of periapsis, for low-inclination orbits we use a related value, the <em>mean longitude</em>, $L=\varpi + M$.</p> <p>The actual angle between the orbiting body and periapsis is called the <em>true anomaly</em>, $\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 <em>eccentric anomaly</em> $E$:</p> <p>$$M = E - e \sin E$$</p> <p>This is called <a href="https://en.wikipedia.org/wiki/Kepler%27s_equation"><strong>Kepler's equation</strong></a>, and it cannot be solved analytically. Once we have $E$ though, there is a relatively simple expression for $\nu$.</p> <h2>Computing Position from Orbital Elements</h2> <p>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.</p> <p>I'll assume that you're using a set of elements like <a href="http://ssd.jpl.nasa.gov/?planet_pos">these from JPL's web site</a>. 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. </p> <p>Calculate the time $t$ in centuries from J2000:</p> <pre><code>// 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); </code></pre> <p>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&ndash;2500):</p> <pre><code>// a0 = 1.00000261; adot = 0.00000562 var a = a0 + adot * tCenturiesFromJ2000; </code></pre> <p>(Note that the values are actually given for "EM <a href="https://en.wikipedia.org/wiki/Barycenter">Barycenter</a>," 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.)</p> <p>Table 2a gives elements that are accurate from 3000 BC to 3000 AD; however, if you use the elements from table 2a, you <em>must</em> supplement them with corrections to $L$ from Table 2b! For example, here is computing the longitude of Saturn:</p> <pre><code>// 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); </code></pre> <p>We don't need to explicitly compute the mean motion and add it to $L$, since both tables include it in $\dot L$.</p> <p>Now we're ready to compute $M$ and $\omega$ (<code>w</code>):</p> <pre><code>var M = L - p \\ p is the longitude of periapsis var w = p - W \\ W is the longitude of the ascending node </code></pre> <p>On to step 2: we need to solve the Kepler equation:</p> <p>$$M = E - e \sin E$$</p> <p>We can solve this numerically using <a href="https://en.wikipedia.org/wiki/Newton%27s_method"><strong>Newton's method</strong></a>. 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:</p> <p>$$E_{i+1} = E_i - f(E_i) / f'(E_i) \\ f'(E) = 1 - e \cos E$$</p> <p>Since the nonlinear part $e \sin E$ is very small, we can start with the estimate $E=M$. Our code looks something like this:</p> <pre><code>E = M; while(true) { var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E)); E -= dE; if( Math.abs(dE) &lt; 1e-6 ) break; } </code></pre> <p>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$:</p> <pre><code>var P = a * (Math.cos(E) - e); var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2)); </code></pre> <p>(<code>P</code> and <code>Q</code> form a 2d coordinate system in the plane of the orbit, with <code>+P</code> pointing towards periapsis.)</p> <p>Finally, we can rotate these coordinates into the full 3d coordinate system:</p> <pre><code>// 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; </code></pre> <p>(<code>x</code>, <code>y</code>, and <code>z</code> will be in units of AU.)</p> <p>And you're done!</p> <hr> <p>A few tips:</p> <ul> <li><p>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:</p> <pre><code>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)); </code></pre> <p>Note that the velocities will be in AU per century.</p></li> <li><p>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!</p></li> <li>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$.</li> </ul> <hr> <p>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 <em>Fundamentals of Astrodynamics</em> by Bate, Mueller, and White <a href="http://www.fgg.uni-lj.si/~/mkuhar/Zalozba/Fundamentals_of_Astrodynamics-Bate_Mueller&amp;White-1971.pdf">(pdf)</a>. 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, <em>Position and Velocity as a Function of Time.</em></p> https://space.stackexchange.com/questions/8911/-/8917#8917 8 Answer by Mark Adler for Determining orbital position at a future point in time Mark Adler https://space.stackexchange.com/users/265 2015-04-26T17:05:15Z 2015-04-26T17:05:15Z <p>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:</p> <p>$$x(t)=a\cos\left({2\pi \left(t-t_0\right)\over T}\right)$$</p> <p>$$y(t)=a\sin\left({2\pi \left(t-t_0\right)\over T}\right)$$</p> <p>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.</p> <p>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:</p> <p>$$T=2\pi\sqrt{a^3\over\mu}$$</p> https://space.stackexchange.com/questions/8911/-/31431#31431 3 Answer by lamont for Determining orbital position at a future point in time lamont https://space.stackexchange.com/users/14420 2018-10-17T20:22:34Z 2018-10-18T08:24:02Z <p>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.</p> <hr> <p>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 <a href="https://en.wikipedia.org/wiki/Universal_variable_formulation" rel="nofollow noreferrer">universal variables</a> formulations of solving the Kepler problem.</p> <p>The best solutions to that are probably going to be Goodyear's method:</p> <p>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 <a href="https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19660027556.pdf" rel="nofollow noreferrer">NASA NTRS TD document on the same material</a>)</p> <p>Shepperd's method:</p> <p>Shepperd, S.W. Celestial Mechanics (1985) 35: 129. <a href="https://doi.org/10.1007/BF01227666" rel="nofollow noreferrer">https://doi.org/10.1007/BF01227666</a></p> <p>Or Danby-Stumpff:</p> <p>Danby, J.M.A. Celestial Mechanics (1987) 40: 303. <a href="https://doi.org/10.1007/BF01235847" rel="nofollow noreferrer">https://doi.org/10.1007/BF01235847</a></p> <p>There is some <a href="https://www.mathworks.com/matlabcentral/fileexchange/48723-matlab-functions-for-two-body-orbit-propagation?focused=3853951&amp;tab=function" rel="nofollow noreferrer">MATLAB code here</a> 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).</p> <p>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 <a href="https://en.wikipedia.org/wiki/Symplectic_integrator" rel="nofollow noreferrer">Symplectic Integration</a>. 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 <a href="https://en.wikipedia.org/wiki/Runge%27s_phenomenon" rel="nofollow noreferrer">Runge's phenomenon</a>, but I haven't wrestled with that, so don't know which approach to take there).</p> <p>If you're going to use Runge-Kutta then <a href="https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method" rel="nofollow noreferrer">Dormand-Prince</a> with dynamic step side and its 3rd order interpolant will be very convenient, and is what Matlab uses in its ode45 solver.</p> <p>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).</p>