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!
Best Answer
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)) $$