from_state_vector inverts arg_pe when Z coordinate is extremely small and negative #44
Description
Reproduction case:
from orbital import Position, Velocity, KeplerianElements, earth
R = Position(0, 2500000, -0.01)
V = Velocity(-16703.901013, 0, 0)
print(KeplerianElements.from_state_vector(R, V, body=earth))Expected output:
KeplerianElements:
Semimajor axis (a) = 10000.000 km
Eccentricity (e) = 0.750000
Inclination (i) = 0.0 deg
Right ascension of the ascending node (raan) = 0.0 deg
Argument of perigee (arg_pe) = 90.0 deg
Mean anomaly at reference epoch (M0) = 0.0 deg
Period (T) = 2:45:52.014055
Reference epoch (ref_epoch) = J2000.000
Mean anomaly (M) = 0.0 deg
Time (t) = 0:00:00
Epoch (epoch) = J2000.000
Actual output:
KeplerianElements:
Semimajor axis (a) = 10000.000 km
Eccentricity (e) = 0.750000
Inclination (i) = 0.0 deg
Right ascension of the ascending node (raan) = 0.0 deg
Argument of perigee (arg_pe) = 270.0 deg
Mean anomaly at reference epoch (M0) = 0.0 deg
Period (T) = 2:45:52.014055
Reference epoch (ref_epoch) = J2000.000
Mean anomaly (M) = 0.0 deg
Time (t) = 0:00:00
Epoch (epoch) = J2000.000
(Note: arg_pe is 270°, it should be 90°)
Analysis:
The incorrectness of the output can be seen by changing the R.z of -0.01 to 0.0, which then changes arg_pe (correctly) to 90°. Such a small change in the position should not rotate the body 180° around its orbit. Alternatively, changing R.z to 1 (correctly) causes raan to become 180° and arg_pe 270°. Either of these would be fine - it doesn't matter what raan and arg_pe are when the inclination is very small, but they have to be 90° apart.
The problem is this code in utilities.elements_from_state_vector:
if ev.z < 0:
arg_pe = 2 * pi - arg_peThis is executed whenever e is not close to 0, inverting arg_pe to compensate for the fact that raan also flips around when the body goes below the horizontal. That is the correct behaviour when i is also not close to 0, but when i is close to 0, raan is hard-coded to 0, so it is not correct to flip arg_pe. So the above code should only run in the "e is not small and i is also not small" case.
Note that if i is exactly 0, this isn't an issue, since ev.z will also be 0. The bug is when the orbit is very very slightly inclined, but i counts as close enough to 0.
I am working on a fix for this alongside several other bugs in from_state_vector and will send a PR soon.