Repository: skyfielders/python-skyfield

Issue:
Getting "ValueError: If eccentricity is 0, argument of periapsis must be 0" in keplerlib.py
Running the following:

```
#!/usr/bin/env python3

import datetime, io, skyfield.api, skyfield.constants, skyfield.data.mpc

print( "Skyfield:", skyfield.__version__ )
latitude = -33
longitude = 151
now = datetime.datetime.strptime( "2020-08-12", "%Y-%m-%d" )
data = "J94T00G  7.0   0.15 J949P   0.00000  353.02318   15.50983    6.76386  0.0000000  0.00358836  42.2543833  E MPC 24084     8   1    3 days              Marsden    0000         1994 TG"
timeScale = skyfield.api.load.timescale( builtin = True )
topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude )
ephemeris = skyfield.api.load( "de421.bsp" )
sun = ephemeris[ "sun" ]
earth = ephemeris[ "earth" ]
with io.BytesIO( data.encode() ) as f:
    dataframe = skyfield.data.mpc.load_mpcorb_dataframe( f ).set_index( "designation", drop = False )
    body = sun + skyfield.data.mpc.mpcorb_orbit( dataframe.loc[ "1994 TG" ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )
```

I get the exception:

```
Traceback (most recent call last):
  File "minorPlanetTest.py", line 17, in <module>
    body = sun + skyfield.data.mpc.mpcorb_orbit( dataframe.loc[ "1994 TG" ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )
  File "/usr/local/lib/python3.6/dist-packages/skyfield/data/mpc.py", line 104, in mpcorb_orbit
    row.designation,
  File "/usr/local/lib/python3.6/dist-packages/skyfield/keplerlib.py", line 173, in _from_mean_anomaly
    gm_au3_d2,
  File "/usr/local/lib/python3.6/dist-packages/skyfield/keplerlib.py", line 310, in ele_to_vec
    raise ValueError('If eccentricity is 0, argument of periapsis must be 0')
ValueError: If eccentricity is 0, argument of periapsis must be 0
```

The minor planet data is sourced from a larger [file](https://minorplanetcenter.net/iau/Ephemerides/Distant/Soft00Distant.txt).

Also note that I am using patched versions of several Skyfield files as a result of other issues.  The files are:


```magnitudelib.py```
```
# -*- coding: utf-8 -*-
"""Routines for computing magnitudes.

Planetary routines adapted from:

https://arxiv.org/pdf/1808.01973.pdf

Which links to:

https://sourceforge.net/projects/planetary-magnitudes/

Which has directories with three successive versions of their magnitude
computation, the most recent of which provides the files on which this
Python code is based:

Ap_Mag_V3.f90
Ap_Mag_Output_V3.txt
Ap_Mag_Input_V3.txt

* ``r`` planet’s distance from the Sun.
* ``delta`` from Earth?
* ``ph_ang`` illumination phase angle (degrees)

"""
from numpy import log10
from .constants import RAD2DEG
from .functions import angle_between, length_of
from .naifcodes import _target_name

def planetary_magnitude(position):
    """Given the position of a planet, return its visual magnitude.

    This prototype function — which so far only supports Mercury, Venus,
    Earth, Jupiter, and Uranus — computes the visual magnitude of a
    planet, given its position relative to an observer.

    >>> from skyfield.api import load
    >>> from skyfield.magnitudelib import planetary_magnitude
    >>> ts = load.timescale()
    >>> t = ts.utc(2020, 7, 31)
    >>> eph = load('de421.bsp')
    >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter'])
    >>> print('%.2f' % planetary_magnitude(astrometric))
    -2.73

    """
    target = position.target
    function = _FUNCTIONS.get(target)
    if function is None:
        name = _target_name(target)
        raise ValueError('cannot compute the magnitude of target %s' % name)

    # Shamelessly treat the Sun as sitting at the Solar System Barycenter.
    sun_to_observer = position.observer_data.bcrs_position
    observer_to_planet = position.position.au
    sun_to_planet = sun_to_observer + observer_to_planet

    r = length_of(sun_to_planet)
    delta = length_of(observer_to_planet)
    ph_ang = angle_between(-sun_to_planet, -observer_to_planet) * RAD2DEG

    return function(r, delta, ph_ang)

def _mercury_magnitude(r, delta, ph_ang):
    distance_mag_factor = 5 * log10(r * delta)
    ph_ang_factor = (
        6.3280e-02 * ph_ang
        - 1.6336e-03 * ph_ang**2
        + 3.3644e-05 * ph_ang**3
        - 3.4265e-07 * ph_ang**4
        + 1.6893e-09 * ph_ang**5
        - 3.0334e-12 * ph_ang**6
    )
    return -0.613 + distance_mag_factor + ph_ang_factor

def _venus_magnitude(r, delta, ph_ang):
    distance_mag_factor = 5 * log10(r * delta)
    if ph_ang < 163.7:
        ph_ang_factor = (
            -1.044E-03 * ph_ang
            + 3.687E-04 * ph_ang**2
            - 2.814E-06 * ph_ang**3
            + 8.938E-09 * ph_ang**4
        )
    else:
        ph_ang_factor = (
            236.05828 + 4.384
            - 2.81914E+00 * ph_ang
            + 8.39034E-03 * ph_ang**2
        )
    return -4.384 + distance_mag_factor + ph_ang_factor

def _earth_magnitude(r, delta, ph_ang):
    distance_mag_factor = 5 * log10 (r * delta)
    ph_ang_factor = -1.060e-03 * ph_ang + 2.054e-04 * ph_ang**2
    return -3.99 + distance_mag_factor + ph_ang_factor

def _jupiter_magnitude(r, delta, ph_ang):
    distance_mag_factor = 5 * log10(r * delta)
    geocentric_phase_angle_limit = 12.0

    if ph_ang <= geocentric_phase_angle_limit:  # TODO: time arrays
        ph_ang_factor = -3.7E-04 * ph_ang + 6.16E-04 * ph_ang**2
    else:
        ph_ang_factor = -2.5 * log10(
            1.0 - 1.507 * (ph_ang / 180.)
            - 0.363 * (ph_ang / 180.)**2
            - 0.062 * (ph_ang / 180.)**3
            + 2.809 * (ph_ang / 180.)**4
            - 1.876 * (ph_ang / 180.)**5
        )

    if ph_ang <= geocentric_phase_angle_limit:
        ap_mag = -9.395 + distance_mag_factor + ph_ang_factor
    else:
        ap_mag = -9.428 + distance_mag_factor + ph_ang_factor

    return ap_mag

def _uranus_magnitude(r, delta, ph_ang,
                      sun_sub_lat_planetog=0.0, earth_sub_lat_planetog=0.0):
    distance_mag_factor = 5.0 * log10 (r * delta)
    sub_lat_planetog = (abs(sun_sub_lat_planetog)
                        + abs(earth_sub_lat_planetog)) / 2.0
    sub_lat_factor = -0.00084 * sub_lat_planetog
    geocentric_phase_angle_limit = 3.1
    ap_mag = -7.110 + distance_mag_factor + sub_lat_factor
    if ph_ang > geocentric_phase_angle_limit:  # TODO: time arrays
        ap_mag += 6.587e-3 * ph_ang + 1.045e-4 * ph_ang**2
    return ap_mag

_FUNCTIONS = {
    199: _mercury_magnitude,
    299: _venus_magnitude,
    399: _earth_magnitude,
    599: _jupiter_magnitude,
    799: _uranus_magnitude,

    # Some planets can be reasonably identified with their barycenter.
    1: _mercury_magnitude,
    2: _venus_magnitude,
    5: _jupiter_magnitude,
    7: _uranus_magnitude,
}
```


```keplerlib.py```
```
from __future__ import division

import sys
import math
from numpy import(abs, amax, amin, arange, arccos, arctan, array, cos, cosh,
                  cross, exp, log, ndarray, newaxis, ones_like, pi, power,
                  repeat, sin, sinh, sqrt, sum, tan, tanh, tile, zeros_like)

from skyfield.constants import AU_KM, DAY_S, DEG2RAD
from skyfield.functions import dots, length_of, mxv
from skyfield.descriptorlib import reify
from skyfield.elementslib import OsculatingElements, normpi
from skyfield.units import Distance, Velocity
from skyfield.vectorlib import VectorFunction

_CONVERT_GM = DAY_S * DAY_S / AU_KM / AU_KM / AU_KM

class _KeplerOrbit(VectorFunction):
    def __init__(self,
                 position,
                 velocity,
                 epoch,
                 mu_au3_d2,
                 center=None,
                 target=None,
        ):
        """ Calculates the position of an object using 2 body propagation

        Parameters
        ----------
        position : Distance
            Position vector at epoch with shape (3,)
        velocity : Velocity
            Velocity vector at epoch with shape (3,)
        epoch : Time
            Time corresponding to `position` and `velocity`
        mu_au_d : float
            Value of mu (G * M) in au^3/d^2
        center : int
            NAIF ID of the primary body, 399 for geocentric orbits, 10 for
            heliocentric orbits
        target : int
            NAIF ID of the secondary body
        """
        self.position_at_epoch = position
        self.velocity_at_epoch = velocity
        self.epoch = epoch
        self.mu_au3_d2 = mu_au3_d2
        self.center = center
        self.target = target

        self._rotation = None  # TODO: make argument?

    @classmethod
    def _from_true_anomaly(cls, p, e, i, Om, w, v,
                          epoch,
                          mu_km_s=None,
                          mu_au3_d2=None,
                          center=None,
                          target=None,
        ):
        """ Creates a `KeplerOrbit` object from elements using true anomaly

        Parameters
        ----------
        p : Distance
            Semi-Latus Rectum
        e : float
             Eccentricity
        i : Angle
            Inclination
        Om : Angle
            Longitude of Ascending Node
        w : Angle
            Argument of periapsis
        v : Angle
            True anomaly
        epoch : Time
            Time corresponding to `position` and `velocity`
        mu_km_s : float
            Value of mu (G * M) in km^3/s^2
        mu_au3_d2 : float
            Value of mu (G * M) in au^3/d^2
        center : int
            NAIF ID of the primary body, 399 for geocentric orbits, 10 for
            heliocentric orbits
        target : int
            NAIF ID of the secondary body
        """
        if (mu_km_s and mu_au3_d2) or (not mu_km_s and not mu_au3_d2):
            raise ValueError('Either mu_km_s or mu_au3_d2 should be used, but not both')

        if mu_au3_d2:
            mu_km_s = mu_au3_d2 * AU_KM**3 / DAY_S**2

        position, velocity = ele_to_vec(p.km,
                                        e,
                                        i.radians,
                                        Om.radians,
                                        w.radians,
                                        v.radians,
                                        mu_km_s,
        )
        return cls(Distance(km=position),
                   Velocity(km_per_s=velocity),
                   epoch,
                   mu_km_s,
                   center=center,
                   target=target,
        )


    @classmethod
    def _from_mean_anomaly(
            cls,
            semilatus_rectum_au,
            eccentricity,
            inclination_degrees,
            longitude_of_ascending_node_degrees,
            argument_of_perihelion_degrees,
            mean_anomaly_degrees,
            epoch,
            gm_km3_s2,
            center=None,
            target=None,
        ):
        """ Creates a `KeplerOrbit` object from elements using mean anomaly

        Parameters
        ----------
        p : Distance
            Semi-Latus Rectum
        e : float
             Eccentricity
        i : Angle
            Inclination
        Om : Angle
            Longitude of Ascending Node
        w : Angle
            Argument of periapsis
        M : Angle
            Mean anomaly
        epoch : Time
            Time corresponding to `position` and `velocity`
        mu_km_s : float
            Value of mu (G * M) in km^3/s^2
        mu_au3_d2 : float
            Value of mu (G * M) in au^3/d^2
        center : int
            NAIF ID of the primary body, 399 for geocentric orbits, 10 for
            heliocentric orbits
        target : int
            NAIF ID of the secondary body
        """
        M = DEG2RAD * mean_anomaly_degrees
        gm_au3_d2 = gm_km3_s2 * _CONVERT_GM
        if eccentricity < 1.0:
            E = eccentric_anomaly(eccentricity, M)
            v = true_anomaly_closed(eccentricity, E)
        elif eccentricity > 1.0:
            E = eccentric_anomaly(eccentricity, M)
            v = true_anomaly_hyperbolic(eccentricity, E)
        else: 
            v = true_anomaly_parabolic(semilatus_rectum_au, gm_au3_d2, M)

        pos, vel = ele_to_vec(
            semilatus_rectum_au,
            eccentricity,
            DEG2RAD * inclination_degrees,
            DEG2RAD * longitude_of_ascending_node_degrees,
            DEG2RAD * argument_of_perihelion_degrees,
            v,
            gm_au3_d2,
        )
        return cls(
            Distance(pos),
            Velocity(vel),
            epoch,
            gm_au3_d2,
            center,
            target,
        )

    def _at(self, time):
        """Propagate the KeplerOrbit to the given Time object

        The Time object can contain one time, or an array of times
        """
        pos, vel = propagate(
            self.position_at_epoch.au,
            self.velocity_at_epoch.au_per_d,
            self.epoch.tt,
            time.tt,
            self.mu_au3_d2,
        )
        if self._rotation is not None:
            pos = mxv(self._rotation, pos)
            vel = mxv(self._rotation, vel)
        return pos, vel, None, None


    @reify
    def elements_at_epoch(self):
        return OsculatingElements(self.position_at_epoch,
                                  self.velocity_at_epoch,
                                  self.epoch,
                                  self._mu_km_s,
        )


    def __str__(self):
        ele = self.elements_at_epoch
        if self.target_name:
            return 'KeplerOrbit {0} {1} -> {2} {3}'.format(self.center,
                                                           self.center_name,
                                                           self.target,
                                                           self.target_name,
            )
        else:
            ele = self.elements_at_epoch
            string = 'KeplerOrbit {0} {1} -> q={2:.2}au e={3:.3f} i={4:.1f} Om={5:.1f} w={6:.1f}'
            return string.format(self.center,
                                 self.center_name,
                                 ele.periapsis_distance.au,
                                 ele.eccentricity,
                                 ele.inclination.degrees,
                                 ele.longitude_of_ascending_node.degrees,
                                 ele.argument_of_periapsis.degrees,
            )


    def __repr__(self):
        return '<{0}>'.format(str(self))


def eccentric_anomaly(e, M):
    """ Iterates to solve Kepler's equation to find eccentric anomaly

    Based on the algorithm in section 8.10.2 of the Explanatory Supplement
    to the Astronomical Almanac, 3rd ed.
    """
    M = normpi(M)
    E = M + e*sin(M)

    max_iters = 100
    iters = 0
    while iters < max_iters:
        dM = M - (E - e*sin(E))
        dE = dM/(1 - e*cos(E))
        E = E + dE
        if abs(dE) < 1e-14: return E
        iters += 1
    else:
        raise ValueError('Failed to converge')


def true_anomaly_hyperbolic(e, E):
    """Calculates true anomaly from eccentricity and eccentric anomaly.

    Equations from the relevant Wikipedia entries.

    """
    return 2.0 * arctan(sqrt((e + 1.0) / (e - 1.0)) * tanh(E/2))


def true_anomaly_closed(e, E):
    """Calculates true anomaly from eccentricity and eccentric anomaly.

    Equations from the relevant Wikipedia entries.

    """
    return 2.0 * arctan(sqrt((1.0 + e) / (1.0 - e)) * tan(E/2))


def true_anomaly_parabolic(p, gm, M):
    """Calculates true anomaly for parabolic orbits from semi-latus rectum, gm, and mean anomaly.

    Equations from https://en.wikipedia.org/wiki/Parabolic_trajectory.

    """
    delta_t = sqrt(2 * p**3 / gm) * M # from http://www.bogan.ca/orbits/kepler/orbteqtn.html
    periapsis_distance = p / 2
    A = 3 / 2 * sqrt(gm / (2 * periapsis_distance**3)) * delta_t
    B = (A + (A**2 + 1))**(1/3)
    return 2 * arctan(B - 1/B)


def ele_to_vec(p, e, i, Om, w, v, mu):
    """Calculates state vectors from orbital elements. Also checks for invalid
    sets of elements.

    Based on equations from this document:

    https://web.archive.org/web/*/http://ccar.colorado.edu/asen5070/handouts/kep2cart_2002.doc
    """
    # Checks that longitude of ascending node is 0 if inclination is 0
    if isinstance(i, ndarray) or isinstance(Om, ndarray):
        if ((i==0)*(Om!=0)).any():
            raise ValueError('If inclination is 0, longitude of ascending node must be 0')
    else:
        if i==0 and Om!=0:
            raise ValueError('If inclination is 0, longitude of ascending node must be 0')

    # Checks that argument of periapsis is 0  if eccentricity is 0
    if isinstance(e, ndarray) or isinstance(w, ndarray):
        if ((e==0)*(w!=0)).any():
            raise ValueError('If eccentricity is 0, argument of periapsis must be 0')
    else:
        if e==0 and w!=0:
            raise ValueError('If eccentricity is 0, argument of periapsis must be 0')

    # Checks that true anomaly is less than arccos(-1/e) for hyperbolic orbits
    if isinstance(e, ndarray) and isinstance(v, ndarray):
        inds = (e>1)
        if (v[inds]>arccos(-1/e[inds])).any():
            raise ValueError('If eccentricity is >1, abs(true anomaly) cannot be more than arccos(-1/e)')
    elif isinstance(e, ndarray) and not isinstance(v, ndarray):
        inds = (e>1)
        if (v>arccos(-1/e[inds])).any():
            raise ValueError('If eccentricity is >1, abs(true anomaly) cannot be more than arccos(-1/e)')
    elif isinstance(v, ndarray) and not isinstance(e, ndarray):
        if e>1 and (v>arccos(-1/e)).any():
            raise ValueError('If eccentricity is >1, abs(true anomaly) cannot be more than arccos(-1/e)')
    else:
        if e>1 and v>arccos(-1/e):
            raise ValueError('If eccentricity is >1, abs(true anomaly) cannot be more than arccos(-1/e)')

    # Checks that inclination is between 0 and pi
    if isinstance(i, ndarray):
        assert ((i>=0) * (i < pi)).all()
    else:
        assert i>=0 and i<pi

    r = p/(1 + e*cos(v))
    h = sqrt(p*mu)
    u = v+w

    X = r*(cos(Om)*cos(u) - sin(Om)*sin(u)*cos(i))
    Y = r*(sin(Om)*cos(u) + cos(Om)*sin(u)*cos(i))
    Z = r*(sin(i)*sin(u))

    X_dot = X*h*e/(r*p)*sin(v) - h/r*(cos(Om)*sin(u) + sin(Om)*cos(u)*cos(i))
    Y_dot = Y*h*e/(r*p)*sin(v) - h/r*(sin(Om)*sin(u) - cos(Om)*cos(u)*cos(i))
    Z_dot = Z*h*e/(r*p)*sin(v) + h/r*sin(i)*cos(u)

    # z and z_dot are independent of Om, so if Om is an array and the other
    # elements are scalars, z and z_dot need to be repeated
    if Z.size!=X.size:
        Z = repeat(Z, X.size)
        Z_dot = repeat(Z_dot, X.size)

    return array([X, Y, Z]), array([X_dot, Y_dot, Z_dot])


dpmax = sys.float_info.max

def bracket(num, end1, end2):
    num[num<end1] = end1
    num[num>end2] = end2
    return num


def find_trunc():
    denom = 2
    factr = 2
    trunc = 1
    x = 1 / denom
    while 1+x > 1:
        denom = denom * (2+factr) * (1+factr)
        factr = factr + 2
        trunc = trunc + 1
        x = 1 / denom
    return trunc


trunc = find_trunc()
odd_factorials = array([math.factorial(i) for i in range(3, trunc*2, 2)])
even_factorials = array([math.factorial(i) for i in range(2, trunc*2, 2)])


def stumpff(x):
    """Calculates Stumpff functions

    Based on the function toolkit/src/spicelib/stmp03.f from the SPICE toolkit,
    which can be downloaded from naif.jpl.nasa.gov/naif/toolkit_FORTRAN.html
    """
    if (x < (-(log(2) + log(dpmax))**2)).any():
        raise ValueError('Argument below lower bound')

    z = sqrt(abs(x))

    c0 = zeros_like(x)
    c1 = zeros_like(x)
    c2 = zeros_like(x)
    c3 = zeros_like(x)

    low = x < -1
    c0[low] = cosh(z[low])
    c1[low] = sinh(z[low])/z[low]
    c2[low] = (1 - c0[low])/x[low]
    c3[low] = (1 - c1[low])/x[low]

    high = x > 1
    c0[high] = cos(z[high])
    c1[high] = sin(z[high])/z[high]
    c2[high] = (1 - c0[high])/x[high]
    c3[high] = (1 - c1[high])/x[high]

    mid = ~low * ~high
    n = sum(mid)
    exponents = tile(arange(0, trunc-1), [n, 1])
    odd_denominators = tile(odd_factorials, [n, 1])
    even_denominators = tile(even_factorials, [n, 1])
    numerators = repeat(x[mid][newaxis].T, trunc-1, axis=1)
    c3[mid] = (sum(power(numerators[:, ::2], exponents[:, ::2])/odd_denominators[:, ::2], axis=1)
        - sum(power(numerators[:, 1::2], exponents[:, 1::2])/odd_denominators[:, 1::2], axis=1))
    c2[mid] = (sum(power(numerators[:, ::2], exponents[:, ::2])/even_denominators[:, ::2], axis=1)
        - sum(power(numerators[:, 1::2], exponents[:, 1::2])/even_denominators[:, 1::2], axis=1))
    c1[mid] = 1 - x[mid]*c3[mid]
    c0[mid] = 1 - x[mid]*c2[mid]

    return c0, c1, c2, c3


def propagate(position, velocity, t0, t1, gm):
    """Propagates a position and velocity vector with an array of times.

    Based on the function toolkit/src/spicelib/prop2b.f from the SPICE toolkit,
    which can be downloaded from naif.jpl.nasa.gov/naif/toolkit_FORTRAN.html

    Parameters
    ----------
    position : ndarray
       Position vector with shape (3,)
    velocity : ndarray
        Velocity vector with shape (3,)
    t0 : float
        Time corresponding to `position` and `velocity`
    t1 : float or ndarray
        Time or times to propagate to
    gm : float
        Gravitational parameter in units that match the other arguments
    """
    if gm <= 0:
        raise ValueError("'gm' should be positive")
    if length_of(velocity) == 0:
        raise ValueError('Velocity vector has zero magnitude')
    if length_of(position) == 0:
        raise ValueError('Position vector has zero magnitude')

    r0 = length_of(position)
    rv = dots(position, velocity)

    hvec = cross(position, velocity)
    h2 = dots(hvec, hvec)

    if h2 == 0:
        raise ValueError('Motion is not conical')

    eqvec = cross(velocity, hvec)/gm + -position/r0
    e = length_of(eqvec)
    q = h2 / (gm * (1+e))

    f = 1 - e
    b = sqrt(q/gm)

    br0 = b * r0
    b2rv = b**2 * rv
    bq = b * q
    qovr0 = q / r0


    maxc = max(abs(br0),
               abs(b2rv),
               abs(bq),
               abs(qovr0/(bq)))

    if f < 0:
        fixed = log(dpmax/2) - log(maxc)
        rootf = sqrt(-f)
        logf = log(-f)
        bound = min(fixed/rootf, (fixed + 1.5*logf)/rootf)
    else:
        logbound = (log(1.5) + log(dpmax) - log(maxc)) / 3
        bound = exp(logbound)

    def kepler(x):
        c0, c1, c2, c3 = stumpff(f*x*x)
        return x*(br0*c1 + x*(b2rv*c2 + x*(bq*c3)))

    dt = t1 - t0

    if not isinstance(dt, ndarray):
        dt = array([dt])
        return_1d_array = True
    else:
        return_1d_array = False

    x = bracket(dt/bq, -bound, bound)
    kfun = kepler(x)

    past = dt < 0
    future = dt > 0

    upper = zeros_like(dt, dtype='float64')
    lower = zeros_like(dt, dtype='float64')
    oldx = zeros_like(dt, dtype='float64')

    lower[past] = x[past]
    upper[future] = x[future]

    while (kfun[past] > dt[past]).any():
        upper[past] = lower[past]
        lower[past] *= 2
        oldx[past] = x[past]
        x[past] = bracket(lower[past], -bound, bound)
        if (x[past] == oldx[past]).any():
            raise ValueError('The input delta time (dt) has a value of {0}.'
                             'This is beyond the range of DT for which we '
                             'can reliably propagate states. The limits for '
                             'this GM and initial state are from {1}'
                             'to {2}.'.format(dt, kepler(-bound), kepler(bound)))
        kfun[past] = kepler(x[past])

    while (kfun[future] < dt[future]).any():
        lower[future] = upper[future]
        upper[future] *= 2
        oldx[future] = x[future]
        x[future] = bracket(upper[future], -bound, bound)
        if (x[future] == oldx[future]).any():
            raise ValueError('The input delta time (dt) has a value of {0}.'
                             'This is beyond the range of DT for which we '
                             'can reliably propagate states. The limits for '
                             'this GM and initial state are from {1} '
                             'to {2}.'.format(dt, kepler(-bound), kepler(bound)))
        kfun[future] = kepler(x[future])

    x = amin(array([upper, amax(array([lower, (lower+upper)/2]), axis=0)]), axis=0)

    lcount = zeros_like(dt)
    mostc = ones_like(dt)*1000
    not_done = (lower < x) * (x < upper)

    while not_done.any():
        kfun[not_done] = kepler(x[not_done])

        high = (kfun > dt) * not_done
        low = (kfun < dt) * not_done
        same = (~high * ~low) * not_done

        upper[high] = x[high]
        lower[low] = x[low]
        upper[same] = lower[same] = x[same]

        condition = not_done * (mostc > 64) * (upper != 0) * (lower != 0)
        mostc[condition] = 64
        lcount[condition] = 0

        # vectorized version of min(upper, max(lower, (upper + lower)/2))
        x[not_done] = amin(array([upper[not_done], amax(array([lower[not_done], (lower[not_done]+upper[not_done])/2]), axis=0)]), axis=0)
        lcount += 1
        not_done = (lower < x) * (x < upper) * (lcount < mostc)

    c0, c1, c2, c3 = stumpff(f*x*x)
    br = br0*c0 + x*(b2rv*c1 + x*(bq*c2))

    pc = 1 - qovr0 * x**2 * c2
    vc = dt - bq * x**3 * c3
    pcdot = -qovr0 / br * x * c1
    vcdot = 1 - bq / br * x**2 * c2

    if return_1d_array:
        position_prop = pc*position + vc*velocity
        velocity_prop = pcdot*position + vcdot*velocity
    else:
        position_prop = pc*tile(position[newaxis].T, dt.size) + vc*tile(velocity[newaxis].T, dt.size)
        velocity_prop = pcdot*tile(position[newaxis].T, dt.size) + vcdot*tile(velocity[newaxis].T, dt.size)

    return position_prop, velocity_prop
```


```mpc.py```
```
# -*- coding: utf-8 -*-
"""Routines for interpreting data from the IAU Minor Planet Center."""

import io
import re

import pandas as pd

from ..data.spice import inertial_frames
from ..keplerlib import _KeplerOrbit
from ..timelib import julian_day

MPCORB_URL = 'https://www.minorplanetcenter.net/iau/MPCORB/MPCORB.DAT.gz'

_MPCORB_COLUMNS = [
    ('designation_packed', (0, 7)),
    ('magnitude_H', (8, 13)),
    ('magnitude_G', (14, 19)),
    ('epoch_packed', (20, 25)),
    ('mean_anomaly_degrees', (26, 35)),
    ('argument_of_perihelion_degrees', (37, 46)),
    ('longitude_of_ascending_node_degrees', (48, 57)),
    ('inclination_degrees', (59, 68)),
    ('eccentricity', (70, 79)),
    ('mean_daily_motion_degrees', (80, 91)),
    ('semimajor_axis_au', (92, 103)),
    ('uncertainty', (105, 106)),
    ('reference', (107, 116)),
    ('observations', (117, 122)),
    ('oppositions', (123, 126)),
    ('observation_period', (127, 136)),
    ('rms_residual_arcseconds', (137, 141)),
    ('coarse_perturbers', (142, 145)),
    ('precise_perturbers', (146, 149)),
    ('computer_name', (150, 160)),
    ('hex_flags', (161, 165)),
    ('designation', (166, 194)),
    ('last_observation_date', (194, 202)),
]
_MPCORB_NECESSARY_COLUMNS = {
    'designation_packed', 'epoch_packed', 'mean_anomaly_degrees',
    'argument_of_perihelion_degrees', 'longitude_of_ascending_node_degrees',
    'inclination_degrees', 'eccentricity', 'mean_daily_motion_degrees',
    'semimajor_axis_au',
}
_MPCORB_DTYPES = {
    # These seem to be ignored by read_fwf()?
    'epoch_packed': 'category',
    'uncertainty': 'category',
    'coarse_perturbers': 'category',
    'precise_perturbers': 'category',
    'computer_name': 'category',
}
_MPCORB_CONVERTERS = {
    'designation_packed': str,
    'hex_flags': str,
}

def load_mpcorb_dataframe(fobj):
    """Parse a Minor Planet Center orbits file into a Pandas dataframe.

    See :doc:`kepler-orbits.rst`.  The MPCORB file format is documented
    at: https://minorplanetcenter.net/iau/info/MPOrbitFormat.html

    """
    # See https://github.com/pandas-dev/pandas/issues/18035
    fobj = io.TextIOWrapper(fobj)
    columns = _MPCORB_COLUMNS
    # if not slow:
    #     keepers = _MPCORB_NECESSARY_COLUMNS
    #     columns = [tup for tup in columns if tup[0] in keepers]
    names, colspecs = zip(*columns)
    df = pd.read_fwf(fobj, colspecs, names=names, dtypes=_MPCORB_DTYPES,
                     converters=_MPCORB_CONVERTERS)
    return df

def mpcorb_orbit(row, ts, gm_km3_s2):
    a = row.semimajor_axis_au
    e = row.eccentricity
    p = a * (1.0 - e*e)

    def n(c):
        return ord(c) - (48 if c.isdigit() else 55)

    def d(s):
        year = 100 * n(s[0]) + int(s[1:3])
        month = n(s[3])
        day = n(s[4])
        return julian_day(year, month, day) - 0.5

    epoch_jd = d(row.epoch_packed)
    t_epoch = ts.tt_jd(epoch_jd)

    minor_planet = _KeplerOrbit._from_mean_anomaly(
        p,
        e,
        row.inclination_degrees,
        row.longitude_of_ascending_node_degrees,
        row.argument_of_perihelion_degrees,
        row.mean_anomaly_degrees,
        t_epoch,
        gm_km3_s2,
        10,
        row.designation,
    )
    minor_planet._rotation = inertial_frames['ECLIPJ2000'].T
    return minor_planet

COMET_URL = 'https://www.minorplanetcenter.net/iau/MPCORB/CometEls.txt'

_COMET_COLUMNS = [
    ('number', (0, 4)),
    ('orbit_type', (4, 5)),
    ('designation_packed', (5, 12)),
    ('perihelion_year', (14, 18)),
    ('perihelion_month', (19, 21)),
    ('perihelion_day', (22, 29)),
    ('perihelion_distance_au', (30, 39)),
    ('eccentricity', (41, 49)),
    ('argument_of_perihelion_degrees', (51, 59)),
    ('longitude_of_ascending_node_degrees', (61, 69)),
    ('inclination_degrees', (71, 79)),
    ('perturbed_epoch_year', (81, 85)),
    ('perturbed_epoch_month', (85, 87)),
    ('perturbed_epoch_day', (87, 89)),
    ('magnitude_H', (91, 95)),
    ('magnitude_G', (96, 100)),
    ('designation', (102, 158)),
    ('reference', (159, 168)),
]
_COMET_FAST_COLUMNS = (
    'perihelion_year', 'perihelion_month', 'perihelion_day',
    'perihelion_distance_au', 'eccentricity', 'argument_of_perihelion_degrees',
    'longitude_of_ascending_node_degrees', 'inclination_degrees',
    'magnitude_H', 'magnitude_G',
    'designation',
)

_fast_comet_re = None
_fast_comet_sub = None

def load_comets_dataframe(fobj):
    """Parse a Minor Planet Center comets file into a Pandas dataframe.

    This imports only the fields essential for computing comet orbits.
    See :func:`~skyfield.data.mpc.load_comets_dataframe_slow()` for a
    slower routine that includes every comet data field.

    See :doc:`kepler-orbits.rst`.  The comet file format is documented
    at: https://www.minorplanetcenter.net/iau/info/CometOrbitFormat.html

    """
    global _fast_comet_re, _fast_comet_sub

    text = fobj.read()

    if _fast_comet_re is None:
        # Build a regular expression that will turn the fixed-width
        # comet file into a CSV that Pandas can import efficiently.

        keepers = set(_COMET_FAST_COLUMNS)
        pat = ['^']
        previous_end = None
        for name, (start, end) in _COMET_COLUMNS:
            if previous_end is not None:
                pat.append(' ' * (start - previous_end))
            keep = name in keepers
            if name == 'designation':
                pat.append('(.*?)  .*')
                break
            else:
                if keep:
                    pat.append('(')
                pat.append('.' * (end - start))
                if keep:
                    pat.append(')')
            previous_end = end

        pat = ''.join(pat)
        sub = ','.join('\\{}'.format(i + 1) for i in range(len(keepers)))

        _fast_comet_re = re.compile(pat.encode('ascii'), re.M)
        _fast_comet_sub = sub.encode('ascii')

    text = _fast_comet_re.sub(_fast_comet_sub, text)
    df = pd.read_csv(io.BytesIO(text), header=None, names=_COMET_FAST_COLUMNS)
    return df

def load_comets_dataframe_slow(fobj):
    """Parse a Minor Planet Center comets file into a Pandas dataframe.

    This routine reads in every single field from the comets data file.
    See :func:`~skyfield.data.mpc.load_comets_dataframe()` for a faster
    routine that omits some of the more expensive comet fields.

    See :doc:`kepler-orbits.rst`.  The comet file format is documented
    at: https://www.minorplanetcenter.net/iau/info/CometOrbitFormat.html

    """
    fobj = io.StringIO(fobj.read().decode('ascii'))
    names, colspecs = zip(*_COMET_COLUMNS)
    df = pd.read_fwf(fobj, colspecs, names=names)
    return df

def comet_orbit(row, ts, gm_km3_s2):
    e = row.eccentricity
    if e == 1.0:
        p = row.perihelion_distance_au / 2
    else:
        a = row.perihelion_distance_au / (1.0 - e)
        p = a * (1.0 - e*e)
    t_perihelion = ts.tt(row.perihelion_year, row.perihelion_month,
                         row.perihelion_day)

    comet = _KeplerOrbit._from_mean_anomaly(
        p,
        e,
        row.inclination_degrees,
        row.longitude_of_ascending_node_degrees,
        row.argument_of_perihelion_degrees,
        0.0,
        t_perihelion,
        gm_km3_s2,
        10,
        row['designation'],
    )
    comet._rotation = inertial_frames['ECLIPJ2000'].T
    return comet

def unpack(designation_packed):
    def n(c):
        return ord(c) - (48 if c.isdigit() else 55)
    s = designation_packed
    s1 = s[1]
    if s1 == '/':
        return s
    return '{0[0]}/{1}{0[2]}{0[3]} {0[4]}{2}{3}'.format(
        s, n(s1), int(s[5:7]), s[7].lstrip('0'))
```


Discussion:
brandon-rhodes: And it makes the code shorter and easier to read! I'm definitely in favor. :slightly_smiling_face: Thanks for expanding the range of orbits this can process.

Patch:
diff --git a/skyfield/keplerlib.py b/skyfield/keplerlib.py
index 8aaab8b..fbc604f 100644
--- a/skyfield/keplerlib.py
+++ b/skyfield/keplerlib.py
@@ -327,22 +327,6 @@ def ele_to_vec(p, e, i, Om, w, v, mu):
 
     https://web.archive.org/web/*/http://ccar.colorado.edu/asen5070/handouts/kep2cart_2002.doc
     """
-    # Checks that longitude of ascending node is 0 if inclination is 0
-    if isinstance(i, ndarray) or isinstance(Om, ndarray):
-        if ((i==0)*(Om!=0)).any():
-            raise ValueError('If inclination is 0, longitude of ascending node must be 0')
-    else:
-        if i==0 and Om!=0:
-            raise ValueError('If inclination is 0, longitude of ascending node must be 0')
-
-    # Checks that argument of periapsis is 0  if eccentricity is 0
-    if isinstance(e, ndarray) or isinstance(w, ndarray):
-        if ((e==0)*(w!=0)).any():
-            raise ValueError('If eccentricity is 0, argument of periapsis must be 0')
-    else:
-        if e==0 and w!=0:
-            raise ValueError('If eccentricity is 0, argument of periapsis must be 0')
-
     # Checks that true anomaly is less than arccos(-1/e) for hyperbolic orbits
     if isinstance(e, ndarray) and isinstance(v, ndarray):
         inds = (e>1)
@@ -359,11 +343,13 @@ def ele_to_vec(p, e, i, Om, w, v, mu):
         if e>1 and v>arccos(-1/e):
             raise ValueError('If eccentricity is >1, abs(true anomaly) cannot be more than arccos(-1/e)')
 
-    # Checks that inclination is between 0 and pi
+    # Checks that inclination is in the range [0, pi]
     if isinstance(i, ndarray):
-        assert ((i>=0) * (i < pi)).all()
+        if not ((i>=0) * (i <= pi)).all():
+            raise ValueError('Inclination outside the range [0, pi] radians')
     else:
-        assert i>=0 and i<pi
+        if not 0 <= i <= pi:
+            raise ValueError('Inclination outside the range [0, pi] radians')
 
     r = p/(1 + e*cos(v))
     h = sqrt(p*mu)


Test patch:
diff --git a/skyfield/tests/test_keplerlib.py b/skyfield/tests/test_keplerlib.py
index 8cbc291..432ec8b 100644
--- a/skyfield/tests/test_keplerlib.py
+++ b/skyfield/tests/test_keplerlib.py
@@ -162,11 +162,21 @@ def test_circular_equatorial():
                 p_eps=1e-2, e_eps=1e-8, i_eps=1e-15)
 
 
+def test_circular_retrograde_equatorial():
+    check_orbit(p=300000, e=0, i=pi, Om=0, w=0, v=1,
+                p_eps=1e-2, e_eps=1e-8, i_eps=1e-15)
+
+
 def test_circular_polar():
     check_orbit(p=300000, e=0, i=pi/2, Om=1, w=0, v=1,
                 p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15)
 
 
+def test_circular_non_zero_arg_of_periapsis():
+    check_orbit(p=300000, e=0, i=.5, Om=1, w=.5, v=1,
+                p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15)
+
+
 def test_elliptical():
     check_orbit(p=300000, e=.3, i=1, Om=0, w=4, v=5,
                 p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-7)
@@ -177,6 +187,11 @@ def test_elliptical_equatorial():
                 p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-7)
 
 
+def test_elliptical_retrograde_equatorial():
+    check_orbit(p=300000, e=.3, i=pi, Om=0, w=4, v=5,
+                p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-7)
+
+
 def test_elliptical_polar():
     check_orbit(p=300000, e=.2, i=pi/2, Om=1, w=2, v=3,
                 p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-8)
@@ -192,6 +207,11 @@ def test_parabolic_equatorial():
                 p_eps=1e-5, e_eps=1e-14, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-13)
 
 
+def test_parabolic_retrograde_equatorial():
+    check_orbit(p=300000, e=1, i=pi, Om=0, w=1, v=2,
+                p_eps=1e-5, e_eps=1e-14, i_eps=1e-15, Om_eps=1e-13, w_eps=1e-13)
+
+
 def test_parabolic_polar():
     check_orbit(p=300000, e=1, i=pi/2, Om=1, w=2, v=3,
                 p_eps=1e-5, e_eps=1e-14, i_eps=1e-14, Om_eps=1e-13, w_eps=1e-13)
@@ -207,6 +227,16 @@ def test_hyperbolic_equatorial():
                 p_eps=1e0, e_eps=1e-6, i_eps=1e-15, Om_eps=1e-15, w_eps=1e-6)
 
 
+def test_hyperbolic_retrograde_equatorial():
+    check_orbit(p=300000, e=1.3, i=pi, Om=0, w=1, v=.5,
+                p_eps=1e0, e_eps=1e-6, i_eps=1e-15, Om_eps=1e-9, w_eps=1e-6)
+
+
 def test_hyperbolic_polar():
     check_orbit(p=300000, e=1.3, i=pi/2, Om=1, w=2, v=.5,
                 p_eps=1e0, e_eps=1e-6, i_eps=1e-10, Om_eps=1e-10, w_eps=1e-6)
+
+
+def test_equatorial_non_zero_longitude_of_ascending_node():
+    check_orbit(p=300000, e=.3, i=0, Om=0, w=4, v=5,
+                p_eps=1e-2, e_eps=1e-8, i_eps=1e-15, w_eps=1e-7)