I am working on a satellite simulator for a project/game and I am stuck on this one bit of physics. So far I have a satellite that revolves around earth on a 2D plane following Keplerian motion using Kepler's equation. Everything is fine and the satellite will orbit many times without problem. I can also change the velocity of the satellite to manipulate apo/peri points as well as augment of periapsis.
However, the problem occurs if I change the velocity once the satellite passes the PI/2 mark, the entire orbit gets flipped (augment of periapsis is offset by PI) and the true anomaly resets to 0.
Here is how it is currently implemented:
With the assumption that the satellite starts at the periapsis at t0, I find the eccentric and true anomaly given a time since periapsis $w$. Then I find the radius at that angle to get the final position.
$$t_0 = 0$$
$$E_0 = 0$$
$$M(0) = E_0 – e \sin E_0$$
$$n = \sqrt{ \frac{GM }{ a^3}}$$
$$E = {\rm kepler}(e, M)$$
$$v = 2 \tan^{-1}{(\sqrt{\frac{(1+e)}{(1-e)}}\tan\frac{E}{2})}$$
$$r = a\frac{1-e^2}{1 + e\cos{v}}$$
When I update the velocity of the satellite (flight path angle $\phi$ is recomputed after change in velocity), I recalculate a new true anomaly using $\phi$, I then find the difference between the old true anomaly and the new true anomaly, and add that to the augment of periapsis $w$. After that I recompute time since periapsis $t$ using the new true anomaly and recompute the orbital parameters as well.
$$vel = \sqrt{GM\frac{2}{r} – \frac{1}{a}}$$
$$\phi = \tan^{-1} ( \frac{e\sin{v}}{1 + e \cos{v}})$$
$$v_2 = \tan^{-1}{\frac{r v g m \cos{\phi} * \sin{\phi}}{r v g m \cos^2{\phi} – 1}}$$
$$w \mathrel{{+}{=}} v_2 – v$$
$$E = \cos^{-1}{ \frac{e + \cos{v}}{1 + e \cos{v}} }$$
$$M = E – e \sin{E}$$
$$n = \sqrt{ \frac{GM}{a^3} }$$
$$t = \frac{M}{n}$$
$${\rm calculateOrbit}()$$
I noticed if I were to $v = v + \pi$ if ever $v_2 < 0$ and $v > 0$, then the weird orbit/position flip won't occur until the $v = \pi$ Otherwise I do not know what the problem is. I would really appreciate it if someone can point me the way, as I have been scratching my head over this for a long time now.
Best Answer
A quick note, your equation for the radius $r$ as a function of the true anomaly is missing the semi-major axis $a$. I prefer to use the symbol $\theta$ for the true anomaly instead of $v$, since I will use that for the total velocity. So: $$ r(\theta)=\frac{a(1-e^2)}{1+e\cos\theta} $$ And the way I would calculate my new trajectory parameters after a velocity change is by evaluating the semi-major axis $a$ and eccentricity $e$: $$ a=\frac{\mu r}{2\mu-rv^2} $$ $$ e=\sqrt{1+\frac{r^3\omega^2}{\mu}\left(\frac{rv^2}{\mu}-2\right)} $$ where $\mu$ is the gravitational parameter (you use $GM$) and $\omega$ is the angular velocity, such that the tangential velocity is equal to $v_{\theta}=r\omega$.
From this you can calculate the true anomaly via: $$ \theta=\cos^{-1}\left(\frac{a(1-e^2)-r}{er}\right) $$ However the arc-cosine will return a value between $0$ and $\pi$ and to cover the rest of the arc (form $-\pi$ to $0$) you have to look at the radial velocity $\dot{r}$ or $v_{r}$. When this is negative it means that you have passed your apoapsis and that the true anomaly will be bigger than $\pi$ or negative (depending on which range you want to use for $\theta$), so the true anomaly would simply become: $$ \theta=-\cos^{-1}\left(\frac{a(1-e^2)-r}{er}\right) $$ PS: This might be a useful question since it contains an explicit expression for the the time as a function of your position.