Determine orbital position after change in velocity

Question

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.

Answer 1

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.