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!
 A: 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 :)
A: 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))
$$
