Why are the true variance of these orbits out by ~pi?

Question

I have an object that is orbit around a point mass in a 2D environment with a known speed, radius and zenith. I calculate the following orbital parameters as such and can confirm that these values are correct (where $R_{p}$ is the periapsis, $R_{a}$ is the apoapsis, $e$ is the eccentricity, $a$ is the semimajor axis, $r$ is the radius, $v$ is the velocity and $GM$ is the product of the mass of the point mass and the gravitational constant):

$C = \frac{2GM}{rv^{2}}$

$R_{p} = r \times \frac{-C+\sqrt{C^{2}-4(1-C)(-\sin^{2}(1-C))}}{2(1-C)}$

$R_{a} = r \times \frac{-C-\sqrt{C^{2}-4(1-C)(-\sin^{2}(1-C))}}{2(1-C)}$

$e = \left | \frac{R_{a} - R_{p}}{R_{a} + R_{p}} \right |$

$a = \frac{R_{a} + R_{p}}{2}$

I then go on to attempt to calculate the true anomaly of the orbit as follows (where $z$ is the zenith and $\theta$ is the true anomaly):

$N = \frac{rv^{2}}{GM}$

$\theta = \tan^{-1} \frac{N\sin z \cos z}{(N\sin^{2} z)-1}$

This is calculates a correct value as long as the radius is greater than a value slightly less than the semiminor axis of the orbit. Adding pi to the true anomaly corrects this error, except where the radius is close in value to the length of the semiminor axis.

Why is the true anomaly off by ~pi in this case?

Answer 1

As indicated in the comments, this is probably a question for Computational Science. But I will try to answer your question in such a manner that it may be on-topic enough to stay here.

I think you can solve this by using Python's atan2(y,x) function, rather than defining $z = y/x$ and using atan(z).

Let me address this concern first:

...except where the radius is close in value to the length of the semiminor axis.

To see what happens when $a$ and $r$ are close in value, let me set them equal to each other. Then plug in the expressions for $R_{p,a}$ into the expression for $a$. That gives $2(C-1)r = r \Rightarrow C = 3/2$.

Also note, from $C = \frac{2GM}{rv^2}$ and $N = \frac{rv^2}{GM}$, that $N = 2/C$ is true regardless of the value of $r$. In the case of $r = a$, this becomes $N = 4/3$.

Plug $N = 4/3$ into your expression for $\theta$. This gives $\theta = \tan^{-1}(\frac{4/3\sin{z}\cos{z}}{(4/3\sin^2{z}) - 1})$. Note that if $z = \pi/3 = 60^\circ$, then the denominator is zero. Mathematically, that's okay, because $\tan^{-1}{\infty}$ is perfectly well-defined. But I think you may be inadvertently dividing by a very small number, and that's making your $\tan^{-1}$ give you unpredictable results. If you're running into machine precision issues, then trying to fix things by adding or subtracting $\pi$ won't help.

The other issue I think you're having is related. If $r$ changes by a bit, then that $4/3$ constant in the denominator of your $\tan^{-1}$ changes, which might change the sign of the expression. If it changes sign, $\tan^{-1}$ might be giving you a value for $\theta$ that's not in the quadrant you expect it to be in. If that happens, it will be off by $\pi$.

If I'm correct about both of those, then atan2 should fix them both. You call atan2 with both the numerator and the denominator, so it has the sign information for both, and knows which quadrant you want. And I think it does a better job of handling values of x close to zero (the folks at Comp Sci would definitely be able to give you a better answer about that than I can).