N = E = +1.0
S = W = -1.0 AU_M = 149597870700             # per IAU 2012 Resolution B2
AU_KM = 149597870.700
ASEC360 = 1296000.0
DAY_S = 86400.0
ASEC2RAD = 4.848136811095359935899141e-6
DEG2RAD = 0.017453292519943296
RAD2DEG = 57.295779513082321
C = 299792458.0                            # m/s
ANGVEL = 7.2921150e-5           # radians/s
ERAD = 6378136.6                # meters
IERS_2010_INVERSE_EARTH_FLATTENING = 298.25642
GS = 1.32712440017987e+20
T0 = 2451545.0
B1950 = 2433282.4235
C_AUDAY = C * DAY_S / AU_M 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,
}
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_g', (91, 95)),
    ('magnitude_k', (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_g', 'magnitude_k',
    'designation', 'reference',
)
def load_mpcorb_dataframe(fobj):
def mpcorb_orbit(row, ts, gm_km3_s2):
    def n(c):
    def d(s):
def load_comets_dataframe(fobj):
def load_comets_dataframe_slow(fobj):
def comet_orbit(row, ts, gm_km3_s2):
def _comet_orbits(rows, ts, gm_km3_s2):
def unpack(designation_packed):
    def n(c): 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_periapsis(
            cls,
            semilatus_rectum_au,
            eccentricity,
            inclination_degrees,
            longitude_of_ascending_node_degrees,
            argument_of_perihelion_degrees,
            t_periapsis,
            gm_km3_s2,
            center=None,
            target=None,
        ):
        """Build a `KeplerOrbit` given its parameters and date of periapsis."""
        gm_au3_d2 = gm_km3_s2 * _CONVERT_GM
        pos, vel = ele_to_vec(
            semilatus_rectum_au,
            eccentricity,
            DEG2RAD * inclination_degrees,
            DEG2RAD * longitude_of_ascending_node_degrees,
            DEG2RAD * argument_of_perihelion_degrees,
            0.0,
            gm_au3_d2,
        )
        return cls(
            Distance(pos),
            Velocity(vel),
            t_periapsis,
            gm_au3_d2,
            center,
            target,
        )

    @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,
                                  mu_km_s = self.mu_au3_d2 / _CONVERT_GM,
        )


    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 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
    """
    gm = atleast_1d(gm)
    if (gm <= 0).any():
        raise ValueError("'gm' should be positive")
    if (length_of(velocity)).any() == 0:
        raise ValueError('Velocity vector has zero magnitude')
    if (length_of(position)).any() == 0:
        raise ValueError('Position vector has zero magnitude')

    if position.ndim == 1:
        position = position[:, newaxis]
    if velocity.ndim == 1:
        velocity = velocity[:, newaxis]

    r0 = length_of(position)
    rv = dots(position, velocity)

    hvec = _cross(position, velocity)
    h2 = dots(hvec, hvec)

    if (h2 == 0).any():
        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 * b * rv
    bq = b * q
    qovr0 = q / r0


    maxc = amax(array([abs(br0),
               abs(b2rv),
               abs(bq),
               abs(qovr0/bq)]), axis=0)

    hyperbolic = (f<0)
    bound = zeros_like(f)

    fixed = log(dpmax/2) - log(maxc[hyperbolic])
    rootf = sqrt(-f[hyperbolic])
    logf = log(-f[hyperbolic])
    bound[hyperbolic] = amin(array([fixed/rootf, (fixed + 1.5*logf)/rootf]), axis=0)

    logbound = (log(1.5) + log(dpmax) - log(maxc[~hyperbolic])) / 3
    bound[~hyperbolic] = exp(logbound)

    # each of these arrays has 1 entry per orbit, so its shape is (#orbits, 1)
    f = f[:, newaxis]
    bq = bq[:, newaxis]
    b2rv = b2rv[:, newaxis]
    br0 = br0[:, newaxis]
    qovr0 = qovr0[:, newaxis]
    bound = bound[:, newaxis]

    def kepler(x):
        _, c1, c2, c3 = stumpff(f*x*x)
        return x*(br0*c1 + x*(b2rv*c2 + x*bq*c3))

    def kepler_1d(x, orb_inds):
        _, c1, c2, c3 = stumpff(x*x*repeat(f, orb_inds))
        return x*(c1*repeat(br0, orb_inds) + x*(c2*repeat(b2rv, orb_inds) + x*(c3*repeat(bq, orb_inds))))


    t1 = atleast_1d(t1)
    t0 = atleast_1d(t0)
    if len(t0) == 1:
        t0 = repeat(t0, position.shape[1])

    # shape of 2 dimensional arrays from here on out should be (#orbits, len(t1))
    dt = t1 - t0[:, newaxis]

    x = dt/bq
    copyto(x, -bound, where=(x<-bound))
    copyto(x, bound, where=(x>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')

    copyto(lower, x, where=past)
    copyto(upper, x, where=future)

    while (kfun[past] > dt[past]).any():
        copyto(upper, lower, where=past)
        lower[past] *= 2
        copyto(oldx, x, where=past)
        orb_ind = sum(past, axis=1)
        x[past] = clip(lower[past], repeat(-bound, orb_ind), repeat(bound, orb_ind))
        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_1d(x[past], orb_ind)

    while (kfun[future] < dt[future]).any():
        copyto(lower, upper, where=future)
        upper[future] *= 2
        copyto(oldx, x, where=future)
        orb_ind = sum(future, axis=1)
        x[future] = clip(upper[future], repeat(-bound, orb_ind), repeat(bound, orb_ind))
        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_1d(x[future], orb_ind)

    x = copy(upper)
    copyto(x, (upper+lower)/2, where=(lower<=upper))

    lcount = zeros_like(dt)
    mostc = full_like(dt, 1000)
    not_done = (lower < x) & (x < upper)

    while not_done.any():
        orb_inds = sum(not_done, axis=1)
        kfun[not_done] = kepler_1d(x[not_done], orb_inds)

        high = (kfun > dt) & not_done
        low = (kfun < dt) & not_done
        same = (~high & ~low) & not_done

        copyto(upper, x, where=(high|same))
        copyto(lower, x, where=(low|same))

        condition = not_done & (mostc > 64) & (upper != 0) & (lower != 0)
        mostc[condition] = 64
        lcount[condition] = 0

        copyto(x, upper, where=(not_done & (lower>upper)))
        copyto(x, (upper+lower)/2, where=(not_done & (lower<=upper)))

        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 * x * c2
    vc = dt - bq * x**3 * c3
    pcdot = -qovr0 / br * x * c1
    vcdot = 1 - bq / br * x * x * c2

    position_prop = pc[newaxis, :, :]*position[:, :, newaxis] + vc[newaxis, :, :]*velocity[:, :, newaxis]
    velocity_prop = pcdot[newaxis, :, :]*position[:, :, newaxis] + vcdot[newaxis, :, :]*velocity[:, :, newaxis]

    return squeeze(position_prop), squeeze(velocity_prop) _CONVERT_GM = DAY_S * DAY_S / AU_KM / AU_KM / AU_KM def ts():
def _data_path(filename):
def compare(value, expected_value, epsilon, mod=False):
def check_types(elements, length):
def horizons_dict(elem, units='km_d', ):
def horizons_array(elem, units='km_d', ):
def test_repr(ts):
def test_single_time(ts):
def test_multiple_times(ts):
def test_equatorial_km_d(ts):
def test_equatorial_km_s(ts):
def test_equatorial_au_d(ts):
def test_ecliptic_km_d(ts):
def test_ecliptic_km_s(ts):
def test_ecliptic_au_d(ts):
def test_extreme_ellipse(ts):
def test_slightly_hyperbolic(ts):
def test_periapsis_time(ts):
def check_orbit(p, e, i, Om, w, v, ts):
def test_circular(ts):
def test_circular_equatorial(ts):
def test_circular_polar(ts):
def test_elliptical(ts):
def test_elliptical_equatorial(ts):
def test_elliptical_polar(ts):
def test_parabolic(ts):
def test_parabolic_equatorial(ts):
def test_parabolic_polar(ts):
def test_hyperbolic(ts):
def test_hyperbolic_equatorial(ts):
def test_hyperbolic_polar(ts):
def test_all_types_at_once(ts):
def test_gm_calculation(ts):
ECLIPTIC = inertial_frames['ECLIPJ2000'] class Angle(Unit):

    def __init__(self, angle=None, radians=None, degrees=None, hours=None,
                 preference=None, signed=False):

        if angle is not None:
            if not isinstance(angle, Angle):
                raise ValueError(_instantiation_instructions)
            self.radians = angle.radians
        elif radians is not None:
            self.radians = _to_array(radians)
        elif degrees is not None:
            self._degrees = degrees = _to_array(_unsexagesimalize(degrees))
            self.radians = degrees / 360.0 * tau
        elif hours is not None:
            self._hours = hours = _to_array(_unsexagesimalize(hours))
            self.radians = hours / 24.0 * tau

        self.preference = (preference if preference is not None
                           else 'hours' if hours is not None
                           else 'degrees')
        self.signed = signed

    @classmethod
    def from_degrees(cls, degrees, signed=False):
        degrees = _to_array(_unsexagesimalize(degrees))
        self = cls.__new__(cls)
        self.degrees = degrees
        self.radians = degrees / 360.0 * tau
        self.preference = 'degrees'
        self.signed = signed
        return self

    @reify
    def radians(self):  # Empty property to provide Sphinx docstring.
        """Radians (𝜏 = 2𝜋 in a circle)."""

    @reify
    def _hours(self):
        return self.radians * 24.0 / tau

    @reify
    def _degrees(self):
        return self.radians * 360.0 / tau

    @reify
    def hours(self):
        r"""Hours (24\ |h| in a circle)."""
        if self.preference != 'hours':
            raise WrongUnitError('hours')
        return self._hours

    @reify
    def degrees(self):
        """Degrees (360° in a circle)."""
        if self.preference != 'degrees':
            raise WrongUnitError('degrees')
        return self._degrees

    def arcminutes(self):
        """Return the angle in arcminutes."""
        return self._degrees * 60.0

    def arcseconds(self):
        """Return the angle in arcseconds."""
        return self._degrees * 3600.0

    def mas(self):
        """Return the angle in milliarcseconds."""
        return self._degrees * 3600000.0

    def __str__(self):
        size = self.radians.size
        if size == 0:
            return 'Angle []'
        if self.preference == 'degrees':
            v = self._degrees
            fmt = _dsgn.format if self.signed else _dfmt.format
            places = 1
        else:
            v = self._hours
            fmt = _hfmt.format
            places = 2
        if size >= 2:
            return '{0} values from {1} to {2}'.format(
                len(v), _sfmt(fmt, v[0], places), _sfmt(fmt, v[-1], places))
        return _sfmt(fmt, v, places)

    def __repr__(self):
        if self.radians.size == 0:
            return '<{0} []>'.format(type(self).__name__)
        else:
            return '<{0} {1}>'.format(type(self).__name__, self)

    def hms(self, warn=True):
        """Convert to a tuple (hours, minutes, seconds).

        All three quantities will have the same sign as the angle itself.

        """
        if warn and self.preference != 'hours':
            raise WrongUnitError('hms')
        sign, units, minutes, seconds = _sexagesimalize_to_float(self._hours)
        return sign * units, sign * minutes, sign * seconds

    def signed_hms(self, warn=True):
        """Convert to a tuple (sign, hours, minutes, seconds).

        The ``sign`` will be either +1 or -1, and the other quantities
        will all be positive.

        """
        if warn and self.preference != 'hours':
            raise WrongUnitError('signed_hms')
        return _sexagesimalize_to_float(self._hours)

    def hstr(self, places=2, warn=True, format=_hfmt):
        """Return a string like ``12h 07m 30.00s``; see `Formatting angles`.

        .. versionadded:: 1.39

           Added the ``format=`` parameter.

        """
        if warn and self.preference != 'hours':
            raise WrongUnitError('hstr')
        hours = self._hours
        shape = getattr(hours, 'shape', ())
        fmt = format.format  # `format()` method of `format` string
        if shape:
            return [_sfmt(fmt, h, places) for h in hours]
        return _sfmt(fmt, hours, places)

    def dms(self, warn=True):
        """Convert to a tuple (degrees, minutes, seconds).

        All three quantities will have the same sign as the angle itself.

        """
        if warn and self.preference != 'degrees':
            raise WrongUnitError('dms')
        sign, units, minutes, seconds = _sexagesimalize_to_float(self._degrees)
        return sign * units, sign * minutes, sign * seconds

    def signed_dms(self, warn=True):
        """Convert to a tuple (sign, degrees, minutes, seconds).

        The ``sign`` will be either +1 or -1, and the other quantities
        will all be positive.

        """
        if warn and self.preference != 'degrees':
            raise WrongUnitError('signed_dms')
        return _sexagesimalize_to_float(self._degrees)

    def dstr(self, places=1, warn=True, format=None):
        """Return a string like ``181deg 52' 30.0"``; see `Formatting angles`.

        .. versionadded:: 1.39

           Added the ``format=`` parameter.

        """
        if warn and self.preference != 'degrees':
            raise WrongUnitError('dstr')
        degrees = self._degrees
        signed = self.signed
        if format is None:
            format = _dsgn if signed else _dfmt
        fmt = format.format  # `format()` method of `format` string
        shape = getattr(degrees, 'shape', ())
        if shape:
            return [_sfmt(fmt, d, places) for d in degrees]
        return _sfmt(fmt, degrees, places)

    def to(self, unit):
        """Convert this angle to the given AstroPy unit."""
        from astropy.units import rad
        return (self.radians * rad).to(unit)

        # Or should this do:
        from astropy.coordinates import Angle
        from astropy.units import rad
        return Angle(self.radians, rad).to(unit) class Distance(Unit):
    """A distance, stored internally as au and available in other units.

    You can initialize a ``Distance`` by providing a single float or a
    float array as either an ``au=``, ``km=``, or ``m=`` parameter.

    You can access the magnitude of the distance with its three
    attributes ``.au``, ``.km``, and ``.m``.  By default a distance
    prints itself in astronomical units (au), but you can take control
    of the formatting and choice of units yourself using standard Python
    numeric formatting:

    >>> d = Distance(au=1)
    >>> print(d)
    1.0 au
    >>> print('{:.2f} km'.format(d.km))
    149597870.70 km

    """
    _warned = False

    def __init__(self, au=None, km=None, m=None):
        if au is not None:
            self.au = _to_array(au)
            """Astronomical units."""
        elif km is not None:
            self.km = km = _to_array(km)
            self.au = km / AU_KM
        elif m is not None:
            self.m = m = _to_array(m)
            self.au = m / AU_M
        else:
            raise ValueError('to construct a Distance provide au, km, or m')

    @classmethod
    def from_au(cls, au):
        self = cls.__new__(cls)
        self.au = _to_array(au)
        return self

    @reify
    def au(self):  # Empty property to provide Sphinx docstring.
        """Astronomical units (the Earth-Sun distance of 149,597,870,700 m)."""

    @reify
    def km(self):
        """Kilometers (1,000 meters)."""
        return self.au * AU_KM

    @reify
    def m(self):
        """Meters."""
        return self.au * AU_M

    def __str__(self):
        n = self.au
        return ('{0} au' if getattr(n, 'shape', 0) else '{0:.6} au').format(n)

    def __repr__(self):
        return '<{0} {1}>'.format(type(self).__name__, self)

    def length(self):
        """Compute the length when this is an |xyz| vector.

        The Euclidean vector length of this vector is returned as a new
        :class:`~skyfield.units.Distance` object.

        >>> from skyfield.api import Distance
        >>> d = Distance(au=[1, 1, 0])
        >>> d.length()
        <Distance 1.41421 au>

        """
        return Distance(au=length_of(self.au))

    def light_seconds(self):
        """Return the length of this vector in light seconds."""
        return self.m / C

    def to(self, unit):
        """Convert this distance to the given AstroPy unit."""
        from astropy.units import au
        return (self.au * au).to(unit) class Velocity(Unit):
    """A velocity, stored internally as au/day and available in other units.

    You can initialize a ``Velocity`` by providing a float or float
    array to its ``au_per_d=`` parameter.

    """
    _warned = False

    # TODO: consider reworking this class to return a Rate object.

    def __init__(self, au_per_d=None, km_per_s=None):
        if km_per_s is not None:
            self.km_per_s = km_per_s = _to_array(km_per_s)
            self.au_per_d = km_per_s * DAY_S / AU_KM
        elif au_per_d is not None:
            self.au_per_d = _to_array(au_per_d)
        else:
            raise ValueError('to construct a Velocity provide'
                             ' au_per_d or km_per_s')

    @reify
    def au_per_d(self):  # Empty property to provide Sphinx docstring.
        """Astronomical units per day."""

    @reify
    def km_per_s(self):
        """Kilometers per second."""
        return self.au_per_d * AU_KM / DAY_S

    @reify
    def m_per_s(self):
        """Meters per second."""
        return self.au_per_d * AU_M / DAY_S

    def __str__(self):
        n = self.au_per_d
        fmt = '{0} au/day' if getattr(n, 'shape', 0) else '{0:.6} au/day'
        return fmt.format(n)

    def __repr__(self):
        return '<{0} {1}>'.format(type(self).__name__, self)

    def to(self, unit):
        """Convert this velocity to the given AstroPy unit."""
        from astropy.units import au, d
        return (self.au_per_d * au / d).to(unit)
import os
from numpy import pi, seterr, linspace
from skyfield.api import load
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import mpc
from skyfield.keplerlib import _KeplerOrbit as KeplerOrbit, propagate, _CONVERT_GM
from skyfield.tests.test_elementslib import compare, ele_to_vec
from skyfield.units import Angle, Distance, Velocity
    from io import BytesIO
    from StringIO import StringIO as BytesIO
    # See the following files in the Skyfield repository:
    #
    # horizons/ceres-orbital-elements
    # horizons/ceres-position

    ts = load.timescale()
    t = ts.tdb_jd(2458886.500000000)

    a = 2.768873850275102E+00 # A
    e = 7.705857791518426E-02 # EC
    p_au = a * (1 - e*e)   # Wikipedia

    k = KeplerOrbit._from_mean_anomaly(
        semilatus_rectum_au=p_au,
        eccentricity=e,
        inclination_degrees=2.718528770987308E+01,
        longitude_of_ascending_node_degrees=2.336112629072238E+01,
        argument_of_perihelion_degrees=1.328964361683606E+02,
        mean_anomaly_degrees=1.382501360489816E+02,
        epoch=t,
        gm_km3_s2=GM_SUN,
        center=None,
        target=None,
    )
    r, v = k._at(t)[:2]
    sun_au = [
        -0.004105894975783999, 0.006739680703224941, 0.002956344702049446,
    ]
    horizons_au = [
        1.334875927366032E+00, -2.239607658161781E+00, -1.328895183461897E+00,
    ]
    epsilon = Distance(m=0.001).au
    assert abs(r + sun_au - horizons_au).max() < epsilon

def test_minor_planet():
    text = (b'00001    3.4   0.15 K205V 162.68631   73.73161   80.28698'
            b'   10.58862  0.0775571  0.21406009   2.7676569  0 MPO492748'
            b'  6751 115 1801-2019 0.60 M-v 30h Williams   0000      '
            b'(1) Ceres              20190915\n')

    ts = load.timescale()
    t = ts.utc(2020, 6, 17)
    eph = load('de421.bsp')
    df = mpc.load_mpcorb_dataframe(BytesIO(text))
    row = df.iloc[0]

    assert row.designation_packed == '00001'
    assert row.designation == '(1) Ceres'

    ceres = mpc.mpcorb_orbit(row, ts, GM_SUN)
    ra, dec, distance = eph['earth'].at(t).observe(eph['sun'] + ceres).radec()

    assert ceres.target == '(1) Ceres'
    assert abs(ra.hours - 23.1437) < 0.00005
    assert abs(dec.degrees - -17.323) < 0.0005

def test_comet():
    text = (b'    CJ95O010  1997 03 29.6333  0.916241  0.994928  130.6448'
            b'  283.3593   88.9908  20200224  -2.0  4.0  C/1995 O1 (Hale-Bopp)'
            b'                                    MPC106342\n')

    ts = load.timescale()
    t = ts.utc(2020, 5, 31)
    eph = load('de421.bsp')
    e = eph['earth'].at(t)

    for loader in mpc.load_comets_dataframe, mpc.load_comets_dataframe_slow:
        df = loader(BytesIO(text))
        row = df.iloc[0]
        k = mpc.comet_orbit(row, ts, GM_SUN)
        p = e.observe(eph['sun'] + k)
        ra, dec, distance = p.radec()

        # The file authorities/mpc-hale-bopp in the repository is the
        # source of these angles.  TODO: can we tighten this bound and
        # drive it to fractions of an arcsecond?

        ra_want = Angle(hours=(23, 59, 16.6))
        dec_want = Angle(degrees=(-84, 46, 58))
        assert abs(ra_want.arcseconds() - ra.arcseconds()) < 2.0
        assert abs(dec_want.arcseconds() - dec.arcseconds()) < 0.2
        assert abs(distance.au - 43.266) < 0.0005

        assert k.target == 'C/1995 O1 (Hale-Bopp)'

def test_comet_with_eccentricity_of_exactly_one():
    ts = load.timescale()
    t = ts.utc(2020, 8, 13)
    planets = load('de421.bsp')
    earth, sun = planets['earth'], planets['sun']

    data = (b'    CK15A020  2015 08  1.8353  5.341055  1.000000  208.8369  '
            b'258.5042  109.1696            10.5  4.0  C/2015 A2 (PANSTARRS)'
            b'                                    MPC 93587')

    with BytesIO(data) as f:
        df = mpc.load_comets_dataframe(f)

    df = df[df['designation'] == 'C/2015 A2 (PANSTARRS)']
    comet = mpc.comet_orbit(df.iloc[0], ts, GM_SUN)
    ra, dec, distance = earth.at(t).observe(sun + comet).radec()

    # These are exactly the RA and declination from the Minor Planet
    # Center for this comet!  (The RA seconds returned by Skyfield
    # actually say "46.45s", which would round up to 46.5, but what's a
    # tenth of an arcsecond among friends?)
    assert str(ra).startswith('18h 46m 46.4')
    assert str(dec).startswith("-72deg 05' 33.")

# Test various round-trips through the kepler orbit object.

def _data_path(filename):
    return os.path.join(os.path.dirname(__file__), 'data', filename)

def check_orbit(p, e, i, Om, w, v,
                p_eps=None, e_eps=None, i_eps=None, Om_eps=None, w_eps=None, v_eps=None):
    pos0, vel0 = ele_to_vec(p, e, i, Om, w, v, mu)

    pos1, vel1 = propagate(pos0, vel0, 0, times, mu)


    orbit = KeplerOrbit(Distance(km=pos1), Velocity(km_per_s=vel1), dummy_time, mu_au3_d2=mu*_CONVERT_GM)