Elements from state vectors does not work for Orbits with 180° inclination.

[jj314]

[mgiuca]

I'm years late but I found the same bug. You get NaN if invoking the utility function elements_from_state_vector, but a crash with an AssertionError if using KeplerianElements.from_state_vector.

Here is a repro case:

from numpy import sqrt
from orbital import Position, Velocity, elements_from_state_vector, \
                    KeplerianElements, earth
RADIUS = 10000000.0
R = Position(RADIUS, 0, 0)
V = Velocity(0, -sqrt(earth.mu / RADIUS), 0)
print(elements_from_state_vector(R, V, earth.mu))
print(KeplerianElements.from_state_vector(R, V, body=earth))

(i.e. a circular retrograde orbit - the position is +X and the velocity is -Y with the correct velocity to establish a circular orbit.)

This doesn't require a circular orbit to trigger. It happens when the orbit is flat (no Z axis) and retrograde (inclination is exactly 180°).

Expected output:

OrbitalElements(a=np.float64(10000000.0), e=np.float64(0.0), i=np.float64(3.141592653589793), raan=0, arg_pe=0, f=0.0)
KeplerianElements:
    Semimajor axis (a)                           =  10000.000 km
    Eccentricity (e)                             =      0.000000
    Inclination (i)                              =    180.0 deg
    Right ascension of the ascending node (raan) =      0.0 deg
    Argument of perigee (arg_pe)                 =      0.0 deg
    Mean anomaly at reference epoch (M0)         =      0.0 deg
    Period (T)                                   = 2:45:52.014054
    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:284: RuntimeWarning: invalid value encountered in scalar divide
  raan = acos(n.x / norm(n))
orbital/utilities.py:290: RuntimeWarning: invalid value encountered in scalar divide
  arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev)))
orbital/utilities.py:302: RuntimeWarning: invalid value encountered in scalar divide
  f = acos(dot(n, r) / (norm(n) * norm(r)))
OrbitalElements(a=np.float64(10000000.0), e=np.float64(0.0), i=np.float64(3.141592653589793), raan=np.float64(nan), arg_pe=np.float64(nan), f=np.float64(nan))
Traceback (most recent call last):
  File "test-fsv.py", line 8, in <module>
    print(KeplerianElements.from_state_vector(R, V, body=earth))
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "orbital/elements.py", line 139, in from_state_vector
    assert self.M0 == oldM0
           ^^^^^^^^^^^^^^^^
AssertionError

Analysis

I am working on a fix along with some other bugs I've found in from_state_vector.

The problem is with the line in elements_from_state_vector:

if abs(i) < SMALL_NUMBER:

This only triggers when i is close to 0, but it should also go into the special case for non-inclined orbits when i is close to 180°. (However, it isn't as simple as just updating this check, because arg_pe will be calculated incorrectly in the retrograde case.) I will upload a PR fixing this soon.