from_state_vector sets arg_pe to nan if orbit is circular and inclined

[mgiuca]

Reproduction case:

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

Expected output:

KeplerianElements:
    Semimajor axis (a)                           =  10000.000 km
    Eccentricity (e)                             =      0.000000
    Inclination (i)                              =     45.0 deg
    Right ascension of the ascending node (raan) =     90.0 deg
    Argument of perigee (arg_pe)                 =      0.0 deg
    Mean anomaly at reference epoch (M0)         =     90.0 deg
    Period (T)                                   = 2:45:52.014054
    Reference epoch (ref_epoch)                  = J2000.000
        Mean anomaly (M)                         =     90.0 deg
        Time (t)                                 = 0:00:00
        Epoch (epoch)                            = J2000.000

Actual output:

orbital/utilities.py:290: RuntimeWarning: invalid value encountered in scalar divide
  arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev)))
KeplerianElements:
    Semimajor axis (a)                           =  10000.000 km
    Eccentricity (e)                             =      0.000000
    Inclination (i)                              =     45.0 deg
    Right ascension of the ascending node (raan) =     90.0 deg
    Argument of perigee (arg_pe)                 =      nan deg
    Mean anomaly at reference epoch (M0)         =     90.0 deg
    Period (T)                                   = 2:45:52.014054
    Reference epoch (ref_epoch)                  = J2000.000
        Mean anomaly (M)                         =     90.0 deg
        Time (t)                                 = 0:00:00
        Epoch (epoch)                            = J2000.000

(Note: arg_pe = nan)

Analysis:

This is due to the line in utilities.elements_from_state_vector:

arg_pe = acos(dot(n, ev) / (norm(n) * norm(ev)))

(This is only relevant in the inclined case where i != 0.) Since ev is zero for circular orbits, there is no periapsis, and this results in a division by zero. It needs a special case like in the non-inclined case where arg_pe is set arbitrarily to 0.0.

I am working on a fix for this alongside several other bugs in from_state_vector and will send a PR soon.