AU_KM = 149597870.700 DAY_S = 86400.0 DEG2RAD = 0.017453292519943296 def dots(v, u):
    """Given one or more vectors in `v` and `u`, return their dot products.

    This works whether `v` and `u` each have the shape ``(3,)``, or
    whether they are each whole arrays of corresponding x, y, and z
    coordinates and have shape ``(3, N)``.

    """
    return (v * u).sum(axis=0) def length_of(xyz):
    """Given a 3-element array |xyz|, return its length.

    The three elements can be simple scalars, or the array can be two
    dimensions and offer three whole series of x, y, and z coordinates.

    """
    return sqrt((xyz * xyz).sum(axis=0)) def mxv(M, v):
    """Matrix times vector: multiply an NxN matrix by a vector."""
    return einsum('ij...,j...->i...', M, v) class reify(object):
    """Adapted from Pyramid's `reify()` memoizing decorator."""
    def __init__(self, method):
        self.method = method
        update_wrapper(self, method)

    def __get__(self, instance, objtype=None):
        if instance is None:
            return self
        value = self.method(instance)
        instance.__dict__[self.__name__] = value
        return value class OsculatingElements(object):
    """One or more sets of osculating orbital elements.

    An ``OsculatingElements`` object can be initialized with the
    following parameters:

    position : Distance object
        Position vector with shape (3,) or (3, n)
    velocity : Velocity object
        Velocity vector with shape (3,) or (3, n)
    time: Time object
        The times of the position and velocity vectors
    mu_km_s: float
        Gravitational parameter (G*M) in units of km^3/s^2

    """
    def __init__(self, position, velocity, time, mu_km_s):
        if mu_km_s <= 0:
            raise ValueError('`mu_km_s` (the standard gravitational parameter '
                             'in km^3/s^2) must be positive and non-zero')

        self._pos_vec = position.km
        self._vel_vec = velocity.km_per_s
        self.time = time
        self._mu = mu_km_s
        self._mu_km_d = mu_km_s * _DAY_S_SQUARED

        self._h_vec = cross(self._pos_vec, self._vel_vec, 0, 0).T
        self._e_vec = eccentricity_vector(self._pos_vec, self._vel_vec, self._mu)
        self._n_vec = node_vector(self._h_vec)

    @reify
    def apoapsis_distance(self):
        Q = apoapsis_distance(self.semi_latus_rectum.km, self.eccentricity)
        return Distance(km=Q)

    @reify
    def argument_of_latitude(self):
        w = self.argument_of_periapsis.radians
        v = self.true_anomaly.radians
        u = (w+v)%tau
        return Angle(radians=u)

    @reify
    def argument_of_periapsis(self):
        w = argument_of_periapsis(self._n_vec, self._e_vec, self._pos_vec, self._vel_vec)
        return Angle(radians=w)

    @reify
    def eccentric_anomaly(self):
        E = eccentric_anomaly(self.true_anomaly.radians,
                              self.eccentricity,
                              self.semi_latus_rectum.km)
        return Angle(radians=E)

    @reify
    def eccentricity(self):
        return length_of(self._e_vec)

    @reify
    def inclination(self):
        i = inclination(self._h_vec)
        return Angle(radians=i)

    @reify
    def longitude_of_ascending_node(self):
        Om = longitude_of_ascending_node(self.inclination.radians, self._h_vec)
        return Angle(radians=Om)

    @reify
    def longitude_of_periapsis(self):
        Om = self.longitude_of_ascending_node.radians
        w = self.argument_of_periapsis.radians
        lp = (Om + w)%tau
        return Angle(radians=lp)

    @reify
    def mean_anomaly(self):
        M = mean_anomaly(self.eccentric_anomaly.radians, self.eccentricity)
        return Angle(radians=M)

    @reify
    def mean_longitude(self):
        L = (self.longitude_of_ascending_node.radians
              + self.argument_of_periapsis.radians
              + self.mean_anomaly.radians)%tau
        return Angle(radians=L)

    @reify
    def mean_motion_per_day(self):
        n = mean_motion(self.semi_major_axis.km,
                        self._mu)
        return Angle(radians=n*DAY_S)

    @reify
    def periapsis_distance(self):
        q = periapsis_distance(self.semi_latus_rectum.km, self.eccentricity)
        return Distance(km=q)

    @reify
    def periapsis_time(self):
        M = mean_anomaly(
            self.eccentric_anomaly.radians,
            self.eccentricity,
            shift=False,
        )
        tp = time_since_periapsis(
            M,
            self.mean_motion_per_day.radians,
            self.true_anomaly.radians,
            self.semi_latus_rectum.km,
            self._mu_km_d,
        )
        ts = self.time.ts
        times = self.time.tdb - tp
        return ts.tdb(jd=times)

    @reify
    def period_in_days(self):
        P = period(self.semi_major_axis.km, self._mu)
        return P/DAY_S

    @reify
    def semi_latus_rectum(self):
        p = semi_latus_rectum(self._h_vec, self._mu)
        return Distance(km=p)

    @reify
    def semi_major_axis(self):
        a = semi_major_axis(self.semi_latus_rectum.km, self.eccentricity)
        return Distance(km=a)

    @reify
    def semi_minor_axis(self):
        b = semi_minor_axis(self.semi_latus_rectum.km, self.eccentricity)
        return Distance(km=b)

    @reify
    def true_anomaly(self):
        v = true_anomaly(self._e_vec, self._pos_vec, self._vel_vec, self._n_vec)
        return Angle(radians=v)

    @reify
    def true_longitude(self):
        Om = self.longitude_of_ascending_node.radians
        w = self.argument_of_periapsis.radians
        v = self.true_anomaly.radians
        l = (Om + w + v)%tau
        return Angle(radians=l)

    def __repr__(self):
        return '<Elements {0} sets>'.format(self.time.tt.size) def normpi(num):
    return (num + pi)%tau - pi 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) class VectorFunction(object):
    """Given a time, computes a corresponding position."""

    ephemeris = None

    @reify
    def vector_name(self):
        return type(self).__name__

    @reify
    def center_name(self):
        return _jpl_name(self.center)

    @reify
    def target_name(self):
        return _jpl_name(self.target)

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

    def __str__(self):
        if self.target is self:
            return self.target_name
        return self.arrow_str()

    def arrow_str(self):
        return '{0} {1} -> {2}'.format(
            self.vector_name, self.center_name, self.target_name,
        )

    def __add__(self, other):
        if self.target != other.center:
            if other.target == self.center:
                self, other = other, self
            else:
                raise ValueError(
                    "you can only add two vectors"
                    " if the target where one of the vectors ends"
                    " is the center where the other vector starts"
                )

        self_vfs = getattr(self, 'vector_functions', None) or (self,)
        other_vfs = getattr(other, 'vector_functions', None) or (other,)

        return VectorSum(self.center, other.target, self_vfs + other_vfs)

    def __neg__(self):
        return ReversedVector(self)

    def __sub__(self, other):
        if self.center != other.center:
            raise ValueError(
                "you can only subtract two vectors"
                " if they both start at the same center"
            )

        self_vfs = getattr(self, 'vector_functions', None) or (self,)
        other_vfs = getattr(other, 'vector_functions', None) or (other,)
        other_vfs = tuple(reversed([-vf for vf in other_vfs]))

        return VectorSum(other.target, self.target, other_vfs + self_vfs)

    def at(self, t):
        """At time ``t``, compute the target's position relative to the center.

        If ``t`` is an array of times, then the returned position object
        will specify as many positions as there were times.  The kind of
        position returned depends on the value of the ``center``
        attribute:

        * Solar System Barycenter: :class:`~skyfield.positionlib.Barycentric`
        * Center of the Earth: :class:`~skyfield.positionlib.Geocentric`
        * Difference: :class:`~skyfield.positionlib.Geometric`
        * Anything else: :class:`~skyfield.positionlib.ICRF`

        """
        if not isinstance(t, Time):
            raise ValueError('please provide the at() method with a Time'
                             ' instance as its argument, instead of the'
                             ' value {0!r}'.format(t))
        p, v, gcrs_position, message = self._at(t)
        center = self.center
        position = build_position(p, v, t, center, self.target)
        position._ephemeris = self.ephemeris
        position._observer_gcrs_au = gcrs_position
        position.message = message
        return position

    def _observe_from_bcrs(self, observer):
        if self.center != 0:
            raise ValueError('you can only observe() a body whose vector'
                             ' center is the Solar System Barycenter,'
                             ' but this vector has the center {0}'
                             .format(self.center_name))
        return _correct_for_light_travel_time(observer, self)

    def geometry_of(self, other):
        raise DeprecationError(
"""the geometry_of() method has, alas, been deprecated

This old method has been replaced by an improved interface.  If you just
need your software working again, install Skyfield 0.9.1 for a quick fix:

    pip install skyfield==0.9.1

Or, to update your old code, replace each operation that looks like:

    position = boston.geometry_of(satellite).at(t)

with the vector math that was previously hiding inside the old method:

    position = (satellite - boston).at(t)""")

    def topos(self, latitude=None, longitude=None, latitude_degrees=None,
              longitude_degrees=None, elevation_m=0.0, x=0.0, y=0.0):
        raise DeprecationError(
"""the topos() method has, alas, been deprecated

This old method has been replaced by an improved interface.  If you just
need your software working again, install Skyfield 0.9.1 for a quick fix:

    pip install skyfield==0.9.1

Or, to update your old code, replace each operation that looks like:

    boston = earth.topos(...)

with the vector math that was previously hiding inside the old method:

    from skyfield.api import Topos
    boston = earth + Topos(...)""")

    def satellite(self, text):
        raise DeprecationError(
"""the satellite() method has, alas, been deprecated

This old method has been replaced by an improved interface.  If you just
need your software working again, install Skyfield 0.9.1 for a quick fix:

    pip install skyfield==0.9.1

Or, to update your old code, replace each operation that looks like:

    sat = earth.satellite(tle_text)

with the vector math (and the little bit of text manipulation) that was
previously hiding inside the old method:

    from skyfield.api import EarthSatellite
    line1, line2 = tle_text.splitlines()[-2:]
    sat = earth + EarthSatellite(line1, line2)""") def _cross(a, b):
    # Nearly 4x speedup over numpy cross(). TODO: Maybe move to .functions?
    return a[_cross120] * b[_cross201] - a[_cross201] * b[_cross120]
import sys
import math
from numpy import(abs, amax, amin, arange, arccos, arctan, array, atleast_1d,
                  clip, copy, copyto, cos, cosh, exp, full_like, log, ndarray,
                  newaxis, pi, power, repeat, sin, sinh, squeeze, sqrt, sum,
                  tan, tanh, 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
from skyfield.sgp4lib import _cross
        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)