En primer lugar, la ecuación de su semieje mayor me parece correcta.
En segundo lugar, no existe una ecuación explícita de su posición en función del tiempo. Sin embargo, hay una a la inversa, por lo que una ecuación para tiempo en función de su posición :
$$ 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} $$
Para pequeños pasos en el tiempo se podría utilizar un polinomio, ya que se puede encontrar cualquier orden de derivada temporal de la anomalía verdadera $\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. $$ Así que..: $$ \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) $$
Sin embargo, para grandes pasos en el tiempo este no será un método preciso, pero se podría utilizar para construir un método numérico iterativo calculando el tiempo a partir de $\theta(t_0+\Delta t)$ para determinar el nuevo $\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}) $$ En determinados casos, este método podría ser inestable (cuando $|\Delta t_{n+1}|>|\Delta t_n|$ ), en ese caso podrías dividir tu paso temporal en múltiples pasos, o aumentar el orden del polinomio. Probablemente hay un montón de otros métodos numéricos que son mejores / más rápido, pero yo no soy un experto en esta área, y esto debería funcionar y darle una idea de cómo resolver esto.
En el caso de las órbitas cerradas, también se puede reducir el tamaño del paso temporal mediante una operación modulo modificada: $\Delta t_{new}=\Delta t-round\left[\frac{\Delta t}{T}\right]T$ donde $T$ es igual al periodo de la órbita.
Editar I
Yo mismo estuve trasteando con esto en MATLAB y me di cuenta de que a menudo me encontraba con un problema. Creo que la mayoría de las veces porque el tiempo calculado era complejo porque $\theta_{n+1}$ superado el alcance máximo, por ejemplo una órbita hiperbólica nunca podrá alcanzar $\theta=\pi$ . Este intervalo viene determinado por el $\tan^{-1}$ Así que cuando $e\ge 1$ que $\theta_{max}=2\tan^{-1}\left(\frac{1+e}{\sqrt{e^2-1}}\right)$ . Ahora parece funcionar de forma fiable, a menos que intente calcular algo que se acerque a la precisión de la máquina. Por cierto, no sé en qué lenguaje informático estás programando este simulador, pero es muy probable que no sea capaz de manejar números complejos. Así que tendrías que tener esto en cuenta al realizar los cálculos, ya que cada cálculo debería devolver números reales, pero ciertos sub-resultados serán imaginarios.
Edición II
También busqué otros métodos para resolver esto, como Runge-Kutta, sin embargo no pude evitar el comportamiento inestable (cuando el paso de tiempo es demasiado grande). Sin embargo, un método que siempre sería estable es una especie de algoritmo de búsqueda binaria, en la que comprobar si un ángulo $\theta_{n+1}=\theta_n + \Delta\theta$ daría un paso en el tiempo menor que el paso deseado y, si es así, guarde ese valor. Y después de cada iteración el paso $\Delta\theta$ se dividirá por dos. Para ello tiene que elegir la inicial $\Delta\theta$ correctamente, a saber $\Delta\theta_0=\tan^{-1}\left(\frac{1+e}{e^2-1}\right)$ cuando $e\ge 1$ y $\Delta\theta_0=\pi$ cuando $e<1$ . La tasa de convergencia no es tan buena como la de los métodos anteriores, especialmente cuando $e$ es grande y el paso de tiempo deseado dará lugar a un radio muy alto. Por lo tanto, es posible que desee cambiar los métodos allí, pero no sé cuando el otro método será estable. Otra opción podría ser utilizar el primer método, pero realizar los cálculos relativos al radio, $r$ En cambio $\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)) $$
