# Calculating Orbital Vectors in the Future

For the 2D space simulator that I am writing (please note, it is not at all homework), I require formulas that will give me the location and velocity of a spaceship, relative to the parent celestial body, at a specific time in the future.

There is only one planetary mass that provides a gravitational attraction for the spaceship, and the spaceship's mass is negligible. From the initial $\vec{r}$ and $\vec{v}$ relative to the planet, I already know how to compute many things, such as:

• standard gravitation parameter $\mu = G M$
• eccentricity $e$
• specific angular momentum $h$
• semi-major axis $a = \frac{h^2}{\mu(1 - e^2)}$ (is this even correct?)
• longitude of periapsis $\varpi$
• eccentric anomoly $E$ from true anomaly $\theta$ (an angle the describes an offset from $\varpi$), and vice versa

So using all these values, or possible more (or less), what formulas can I use to compute the position and velocity in $t$ seconds? I have already tried using the following formulas to compute the position in the future:

• mean anomaly $M(E) = E - e sin E$
• mean anomaly at periapsis $M_0 = M(E = \theta = 0)$
• mean anomaly at a certain time $M(t) = M_0 + t \sqrt{\frac{\mu}{a^3}}$

However, those formulas don't work for certain orbits (such as hyperbolic orbits, where $a < 0$). Also, I may have just programmed the formulas incorrectly, but using mean anomaly, my simulator does not correctly determine position in the future.

So, I have already tried several approaches to computing position and velocity in the future, and they didn't work. What are the correct formulas?

Thanks!

• There is a lesson to be learned here - simulations (i.e. differential equation solvers) are very often more easily implemented via timestepping rather than via complicated formulas. You may be interested in an answer I gave to What is the correct way of integrating in astronomy simulations? – user10851 Nov 14 '13 at 18:03
• @ChrisWhite I already run the simulation with timesteps, even with small fixes for extra numerical stability, but I know that it is possible to compute exact position/velocity as a function of time for non-decaying orbits. Hence my question... – feralin Nov 15 '13 at 14:10
• Hm, I'm a bit confused: While implementing a small animation which shows an objects position during the orbit given the starting position and velocity I noticed all my objects are moving in the same direction. Then I realized that M(t) does not contain any information about the direction. What's wrong here? – Gigo Nov 19 '13 at 4:13
• I think this might be caused by the reduction from 3D to 2D. In 3D there are more angles to rotate around so the direction of movement results in different angles there. In 2D I now negate the angular-speed (or what the factor within M(t) is called) when the angular momentum h is negative. Can someone comment on that? – Gigo Nov 19 '13 at 4:29

First of your equation of your semi-major axis seems correct to me.

Secondly there is no explicit equation of your position as a function of time. However there is one the other way around, so an equation for time as a function of your position:

$$t(\theta_0,\theta_1)=\sqrt{\frac{a^3}{\mu}}\left[2\tan^{-1}\left(\frac{\sqrt{1-e^2}\tan{\frac{\theta}{2}}}{1+e}\right)-\frac{e\sqrt{1-e^2}\sin{\theta}}{1+e\cos{\theta}}\right]^{\theta_1}_{\theta_0}$$

For small steps in time you could use a polynomial, since you can find any order time derivative of the true anomaly $\theta$: $$\omega(\theta)=\frac{\delta}{\delta t}\theta=\frac{\delta \theta / \delta \theta}{\delta t / \delta \theta}=\frac{1}{\delta t / \delta \theta}=\sqrt{\frac{\mu}{a^3\left(1-e^2\right)^3}}\left(1+e\cos{\theta}\right)^2,$$ $$\alpha(\theta)=\frac{\delta^2}{\delta t^2}\theta=\frac{\delta}{\delta t}\omega=\frac{\delta \omega / \delta \theta}{\delta t / \delta \theta}=\frac{-2e\mu\sin{\theta}}{a^3(1-e^2)^3}(1+e\cos{\theta})^3,$$ $$ect.$$ So: $$\theta(t_0+\Delta t)=\theta_{t_0}+\Delta t\omega(\theta_{t_0})+\frac{1}{2}\Delta t^2\alpha(\theta_{t_0})+O(\Delta t^3)$$

However for big steps in time this will not be an accurate method, but you could use it to construct an iterative numerical method by calculating the time from $\theta(t_0+\Delta t)$ to determine the new $\Delta t$: $$\theta_{n+1}=\theta_n+\Delta t_n\omega(\theta_n)+\frac{1}{2}\Delta t_n^2\alpha(\theta_n)$$ $$\Delta t_{n+1}=\Delta t_n-t(\theta_n, \theta_{n+1})$$ In certain cases it might be that this method would be unstable (when $|\Delta t_{n+1}|>|\Delta t_n|$) in that case you could split up your time step into multiple steps, or increase the order of the polynomial. There are probably a lot of other numerical methods which are better/faster, but I am not an expert in this area, and this should work and give you in idea of how to solve this.

For closed orbits you might also be able to reduce the size of your time step by taking a modified modulo operation: $\Delta t_{new}=\Delta t-round\left[\frac{\Delta t}{T}\right]T$, where $T$ is equal to the period of the orbit.

Edit I
I was messing with this myself in MATLAB and noticed that I often encountered a problem. I think most of the time because the calculated time was complex because $\theta_{n+1}$ exceeded the maximum range, for example a hyperbolic orbit will never be able to reach $\theta=\pi$. This range is determined by the $\tan^{-1}$, so when $e\ge 1$ than $\theta_{max}=2\tan^{-1}\left(\frac{1+e}{\sqrt{e^2-1}}\right)$. Now it does seems to work reliable, unless I try to calculate something which is getting near the machine accuracy. By the way I do not know in which computer language you are programming this simulator, but there is a big change that it will not be able to handle complex numbers. So you would have to keep this in mind when performing calculations, since each calculation should return real numbers, but certain sub-results will imaginary.

Edit II
I also looked into other methods to solve this, such as Runge–Kutta, however I was not able to prevent unstable behavior (when the time step is to big). However a method which would always be stable is a sort of binary search algorithm, in which check if an angle $\theta_{n+1}=\theta_n + \Delta\theta$ would yield a step in time smaller than you desired step and if so save that value. And after each iteration the step $\Delta\theta$ will be divided by two. For this you do have to choose the initial $\Delta\theta$ correctly, namely $\Delta\theta_0=\tan^{-1}\left(\frac{1+e}{e^2-1}\right)$ when $e\ge 1$ and $\Delta\theta_0=\pi$ when $e<1$. The convergent rate is not as good as that of the previous methods, especially when $e$ is big and the desired time step will result in a very high radius. So you might want to switch methods there, but I do not know when the other method will be stable. Another option might be to use the first method, but to perform the calculations relative to the radius, $r$, instead $\theta$: $$t(r)=\sqrt{\frac{a^3}{\mu}}\left(2\tan^{-1}{\sqrt{\frac{r-a(1-e)}{a(1+e)-r}}}-\sqrt{e^2-\left(\frac{r}{a}-1\right)^2}\right)$$ $$\dot{r}=\sqrt{\frac{\mu}{a}}\frac{\sqrt{e^2a^2-(r-a)^2}}{r}$$ $$\ddot{r}=\mu\left(\frac{a(1-e^2)}{r^3}-\frac{1}{r^2}\right)$$ $$\dddot{r}=\sqrt{\frac{\mu^3}{a}}\frac{\sqrt{e^2a^2-(r-a)^2}}{r^5}(2r-3a(1-e^2))$$

• Thank you for all this work! However, I think you misunderstood what I was asking a bit... I won't need to integrate timesteps when in a non-decaying orbit. When I know that the orbit will never pass through a planet's atmosphere, I can simply keep the current orbit and compute (near-)exact position in the future, to determine where the ship is. So, I think all I really need is your first formula. Let me see if it works. – feralin Nov 15 '13 at 14:24
• @feralin any results yet about whether it worked or not? – fibonatic Nov 17 '13 at 23:43

There's no closed form solution for the eccentric anomaly as a function of time. But there are series solutions that converge rather well. This is actually why Bessel functions were discovered! See here or here. You will have to solve the relationship between the eccentric anomaly and mean anomaly, i.e., $M = E - e \cdot \sin E$ numerically, e.g., using Newton's method or the like.

If you take a look at the following plot of $E$ vs $M$ (at left),

you can see that inverting this function can be done very quickly (and very stably) using Newton's method or bisection, since it is monotonic and single-valued.

So once you have $E$ in terms of $M(t) = M_0 - \omega t$ numerically, you can find $r$ and $\theta$,

$r = a(1-e \cos E)$, and $\theta$ which you said you already have.

If you are just interested in speed, use the Vis-Viva equation $v^2 = GM \left({{ 2 \over{r}} - {1 \over{a}}}\right)$ however, with the position of the particle known in polar coordinates, it is just a matter of differentiation to get it's velocity vector, and substituting the eccentric anomaly back in. The meat of it is on the last page of this link, although it's easy to find in most texts on classical mechanics. I won't repeat it here because I will probably just make a typo and it will have to be repeated anyway.

For anyone who is still interested, here is another treatment of the problem in detail:

Convert Kepler elements to Cartesian coordinates

This is pretty much the canonical way of approaching this problem, and has been for hundreds of years.

This method can also be used in the hyperbolic case (hyperbolic Kepler's equation) or in the parabolic case (Baker's equation), although the equations are somewhat different.

If you're simulating things you might want to consider using a numerics library anyway because then you can start to deal with multibody problems like orbital transfers, and of course, thrust :)

• Can you also comment on how to get the velocity at time t? – Gigo Nov 18 '13 at 4:07
• I made the edit, although I suppose I can do the heavy lifting upon request. – lionelbrits Nov 18 '13 at 10:42
• Thanks for the edit! I can calculate v once r is known using the vis-viva equation, but what is the easiest way to get the velocity including its direction? – Gigo Nov 19 '13 at 2:15
• The whole shebang is done on the last page of [bruce-shapiro.com/pair/ElementConversionRecipes.pdf ] (Recipes for Conversion of Orbital Elements) except the inclination $i=0$. It is also [castle.strw.leidenuniv.nl/ca2013/assets/Assignments_1.pdf ](here) and any search for "convert kepler to cartesian". I don't really know why one would do this calculation any other way, other than for exercise, since this is the way astronomers have done it since before Calculus was invented, but I digress :) – lionelbrits Nov 19 '13 at 11:55