from_state_vector sometimes crashes with f=nan at the periapsis #40
Description
Reproduction case:
from orbital import Position, Velocity, elements_from_state_vector, \
KeplerianElements, earth
R = Position(0, -1767766.952966369, -1767766.952966369)
V = Velocity(16703.901013, 0, 0)
print(elements_from_state_vector(R, V, earth.mu))
print(KeplerianElements.from_state_vector(R, V, body=earth))(Unfortunately, this is a super-test-case that triggers #38 and #39 as well; it's the only one I could find that triggers the f=nan case also.)
Expected output:
OrbitalElements(a=np.float64(10000000.000527445), e=np.float64(0.7500000000131862), i=np.float64(0.7853981633974483), raan=0.0, arg_pe=np.float64(4.71238898038469), f=np.float64(0.0))
KeplerianElements:
Semimajor axis (a) = 10000.000 km
Eccentricity (e) = 0.750000
Inclination (i) = 45.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
Actual output:
orbital/utilities.py:311: RuntimeWarning: invalid value encountered in arccos
f = acos(dot(ev, r) / (norm(ev) * norm(r)))
OrbitalElements(a=np.float64(10000000.000527445), e=np.float64(0.7500000000131862), i=np.float64(0.7853981633974483), raan=np.float64(0.0), arg_pe=np.float64(4.71238898038469), f=np.float64(nan))
Traceback (most recent call last):
File "test-fsv.py", line 7, in <module>
print(KeplerianElements.from_state_vector(R, V, body=earth))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/matt/src-pkg/orbital/orbital/elements.py", line 139, in from_state_vector
assert self.M0 == oldM0
^^^^^^^^^^^^^^^^
AssertionError
Analysis:
The problem happens unpredictably on this line:
f = acos(dot(ev, r) / (norm(ev) * norm(r)))When at the periapsis, ev and r are co-linear. Mathematically, dot(ev, r) should be equal to norm(ev) * norm(r) in this case, so this should result in acos(1) = 0.0. However, due to floating point rounding errors, dot(ev, r) can be slightly larger than norm(ev) * norm(r), resulting in acos(1.0000000000000002), which is nan (and results in the warning "invalid value encountered in arccos").
To fix this, clamp the input to acos to the range [0, 1]. I am working on a PR for this and several other issues with from_state_vector, which I will put up soon.