def interp(*args):
    return np.array(np.interp(*args)) 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 class Splines(object):
    def __init__(self, table):
        table = _to_array(table)
        if len(table.shape) < 2:  # Let caller provide a single row.
            table = table.reshape(table.shape + (1,))
        self.table = table
        self.lower = lower = table[0]
        self.upper = upper = table[1]
        self._width = upper - lower
        self._n = arange(len(lower))
        self.coefficients = table[2:]

    def __call__(self, x):
        i = interp(x, self.lower, self._n)
        i = i.astype(int)
        t = (x - self.lower[i]) / self._width[i]
        coefficients = iter(self.coefficients)
        value = next(coefficients)[i]
        for c in coefficients:
            value *= t
            value += c[i]
        return value

    @reify
    def derivative(self):
        columns = [self.table[0], self.table[1]]
        coefficients = self.table[2:-1]
        for i, c in enumerate(coefficients):
            n = len(coefficients) - i
            columns.append(n * c / self._width)
        return Splines(columns) def build_spline_given_ends(x0, y0, slope0, x1, y1, slope1):
    width = x1 - x0
    slope0 = slope0 * width
    slope1 = slope1 * width
    a0 = y0
    a1 = slope0
    a2 = -2*slope0 - slope1 - 3*y0 + 3*y1
    a3 = slope0 + slope1 + 2*y0 - 2*y1
    return x0, x1, a3, a2, a1, a0 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 def sidereal_time(t):
    """Compute Greenwich Mean Sidereal Time (GMST) in hours at time ``t``."""

    theta = earth_rotation_angle(t.whole, t.ut1_fraction)

    # The equinox method.  See Circular 179, Section 2.6.2.
    # Precession-in-RA terms in mean sidereal time taken from third
    # reference, eq. (42), with coefficients in arcseconds.

    t = (t.whole - T0 + t.tdb_fraction) / 36525.0
    st =        ( 0.014506 +
        (((( -    0.0000000368   * t
             -    0.000029956  ) * t
             -    0.00000044   ) * t
             +    1.3915817    ) * t
             + 4612.156534     ) * t)

    return (st / 54000.0 + theta * 24.0) % 24.0 def earth_rotation_angle(jd_ut1, fraction_ut1=0.0):
    """Return the value of the Earth Rotation Angle (theta) for a UT1 date.

    Uses the expression from the note to IAU Resolution B1.8 of 2000.
    Returns a fraction between 0.0 and 1.0 whole rotations.

    """
    th = 0.7790572732640 + 0.00273781191135448 * (jd_ut1 - T0 + fraction_ut1)
    return (th % 1.0 + jd_ut1 % 1.0 + fraction_ut1) % 1.0 def build_matrix():
    def rotation_at(t):
def build_ecliptic_matrix(t):
    def rotation_at(t):
    def rotation_at(t):
    def _dRdt_times_RT_at(t):
    def rotation_at(t):
    def _dRdt_times_RT_at(t):
    def rotation_at(self, t):
    def __init__(self, doc, matrix):
    def rotation_at(self, t):
class ICRS(object):
class true_equator_and_equinox_of_date(object):
class tirs(object):
class itrs(object):
class ecliptic_frame(object):
class InertialFrame(object):
        R = mxm(rot_z(-t.gast * tau / 24.0), t.M)
            R = mxm(t.polar_motion_matrix(), R) A = A() class A(object):
    """Allow literal NumPy arrays to be spelled ``A[1, 2, 3]``."""
    __getitem__ = array def mxm(M1, M2):
    """Matrix times matrix: multiply two NxN matrices."""
    return einsum('ij...,jk...->ik...', M1, M2) def mxmxm(M1, M2, M3):
    """Matrix times matrix times matrix: multiply 3 NxN matrices together."""
    return einsum('ij...,jk...,kl...->il...', M1, M2, M3) def load_bundled_npy(filename):
    """Load a binary NumPy array file that is bundled with Skyfield."""
    data = get_data('skyfield', 'data/{0}'.format(filename))
    return load(BytesIO(data)) def rot_x(theta):
    c = cos(theta)
    s = sin(theta)
    zero = theta * 0.0
    one = zero + 1.0
    return array(((one, zero, zero), (zero, c, -s), (zero, s, c))) def rot_y(theta):
    c = cos(theta)
    s = sin(theta)
    zero = theta * 0.0
    one = zero + 1.0
    return array(((c, zero, s), (zero, one, zero), (-s, zero, c))) def rot_z(theta):
    c = cos(theta)
    s = sin(theta)
    zero = theta * 0.0
    one = zero + 1.0
    return array(((c, -s, zero), (s, c, zero), (zero, zero, one))) def _to_array(value):
    """Convert plain Python sequences into NumPy arrays.

    This helps Skyfield endpoints convert caller-provided tuples and
    lists into NumPy arrays.  If the ``value`` is not a sequence, then
    it is coerced to a Numpy float object, but not an actual array.

    """
    if hasattr(value, 'shape'):
        return value
    elif hasattr(value, '__len__'):
        return array(value)
    else:
        return float64(value) def build_nutation_matrix(mean_obliquity_radians,
                          true_obliquity_radians,
                          psi_radians):
    """Generate the nutation rotation matrix, given three nutation parameters.

    The input angles can be simple floats.  Or, they can be arrays of
    the same length, in which case the output matrix will have an extra
    dimension of that same length providing *n* rotation matrices.

    """
    cobm = cos(mean_obliquity_radians)
    sobm = sin(mean_obliquity_radians)
    cobt = cos(true_obliquity_radians)
    sobt = sin(true_obliquity_radians)
    cpsi = cos(psi_radians)
    spsi = sin(psi_radians)

    return array(((cpsi,
                  -spsi * cobm,
                  -spsi * sobm),
                  (spsi * cobt,
                   cpsi * cobm * cobt + sobm * sobt,
                   cpsi * sobm * cobt - cobm * sobt),
                  (spsi * sobt,
                   cpsi * cobm * sobt - sobm * cobt,
                   cpsi * sobm * sobt + cobm * cobt))) def equation_of_the_equinoxes_complimentary_terms(jd_tt):
    """Compute the complementary terms of the equation of the equinoxes.

    `jd_tt` - Terrestrial Time: Julian date float, or NumPy array of floats

    """
    # Interval between fundamental epoch J2000.0 and current date.

    t = (jd_tt - T0) / 36525.0

    # Build array for intermediate results.

    shape = getattr(jd_tt, 'shape', ())
    fa = zeros((14,) if shape == () else (14, shape[0]))

    # Mean Anomaly of the Moon.

    fa[0] = ((485868.249036 +
              (715923.2178 +
              (    31.8792 +
              (     0.051635 +
              (    -0.00024470)
              * t) * t) * t) * t) * ASEC2RAD
              + (1325.0*t % 1.0) * tau)

    # Mean Anomaly of the Sun.

    fa[1] = ((1287104.793048 +
              (1292581.0481 +
              (     -0.5532 +
              (     +0.000136 +
              (     -0.00001149)
              * t) * t) * t) * t) * ASEC2RAD
              + (99.0*t % 1.0) * tau)

    # Mean Longitude of the Moon minus Mean Longitude of the Ascending
    # Node of the Moon.

    fa[2] = (( 335779.526232 +
              ( 295262.8478 +
              (    -12.7512 +
              (     -0.001037 +
              (      0.00000417)
              * t) * t) * t) * t) * ASEC2RAD
              + (1342.0*t % 1.0) * tau)

    # Mean Elongation of the Moon from the Sun.

    fa[3] = ((1072260.703692 +
              (1105601.2090 +
              (     -6.3706 +
              (      0.006593 +
              (     -0.00003169)
              * t) * t) * t) * t) * ASEC2RAD
              + (1236.0*t % 1.0) * tau)

    # Mean Longitude of the Ascending Node of the Moon.

    fa[4] = (( 450160.398036 +
              (-482890.5431 +
              (      7.4722 +
              (      0.007702 +
              (     -0.00005939)
              * t) * t) * t) * t) * ASEC2RAD
              + (-5.0*t % 1.0) * tau)

    fa[ 5] = (4.402608842 + 2608.7903141574 * t)
    fa[ 6] = (3.176146697 + 1021.3285546211 * t)
    fa[ 7] = (1.753470314 +  628.3075849991 * t)
    fa[ 8] = (6.203480913 +  334.0612426700 * t)
    fa[ 9] = (0.599546497 +   52.9690962641 * t)
    fa[10] = (0.874016757 +   21.3299104960 * t)
    fa[11] = (5.481293872 +    7.4781598567 * t)
    fa[12] = (5.311886287 +    3.8133035638 * t)
    fa[13] = (0.024381750 +    0.00000538691 * t) * t

    fa %= tau

    # Evaluate the complementary terms.

    a = ke1.dot(fa)
    c_terms = se1_0 * sin(a)
    c_terms += se1_1 * cos(a)
    c_terms *= t

    a = ke0_t.dot(fa)
    c_terms += se0_t_0.dot(sin(a))
    c_terms += se0_t_1.dot(cos(a))

    c_terms *= ASEC2RAD
    return c_terms def iau2000a_radians(t, fundamental_argument_terms=5, lunisolar_terms=687,
                     planetary_terms=687):
    """Return the IAU 2000A angles delta-psi and delta-epsilon in radians."""
    d_psi, d_eps = iau2000a(t.tt, fundamental_argument_terms, lunisolar_terms,
                            planetary_terms)
    d_psi *= _TENTH_USEC_2_RAD
    d_eps *= _TENTH_USEC_2_RAD
    return d_psi, d_eps def mean_obliquity(jd_tdb):
    """Return the mean obliquity of the ecliptic in arcseconds.

    `jd_tt` - TDB time as a Julian date float, or NumPy array of floats

    """
    # Compute time in Julian centuries from epoch J2000.0.

    t = (jd_tdb - T0) / 36525.0

    # Compute the mean obliquity in arcseconds.  Use expression from the
    # reference's eq. (39) with obliquity at J2000.0 taken from eq. (37)
    # or Table 8.

    epsilon = (((( -  0.0000000434   * t
                   -  0.000000576  ) * t
                   +  0.00200340   ) * t
                   -  0.0001831    ) * t
                   - 46.836769     ) * t + 84381.406

    return epsilon def compute_precession(jd_tdb):
    """Return the rotation matrices for precessing to an array of epochs.

    `jd_tdb` - array of TDB Julian dates

    The array returned has the shape `(3, 3, n)` where `n` is the number
    of dates that have been provided as input.

    """
    eps0 = 84381.406

    # 't' is time in TDB centuries.

    t = (jd_tdb - T0) / 36525.0

    # Numerical coefficients of psi_a, omega_a, and chi_a, along with
    # epsilon_0, the obliquity at J2000.0, are 4-angle formulation from
    # Capitaine et al. (2003), eqs. (4), (37), & (39).

    psia   = ((((-    0.0000000951  * t
                 +    0.000132851 ) * t
                 -    0.00114045  ) * t
                 -    1.0790069   ) * t
                 + 5038.481507    ) * t

    omegaa = ((((+    0.0000003337  * t
                 -    0.000000467 ) * t
                 -    0.00772503  ) * t
                 +    0.0512623   ) * t
                 -    0.025754    ) * t + eps0

    chia   = ((((-    0.0000000560  * t
                 +    0.000170663 ) * t
                 -    0.00121197  ) * t
                 -    2.3814292   ) * t
                 +   10.556403    ) * t

    eps0 = eps0 * ASEC2RAD
    psia = psia * ASEC2RAD
    omegaa = omegaa * ASEC2RAD
    chia = chia * ASEC2RAD

    sa = sin(eps0)
    ca = cos(eps0)
    sb = sin(-psia)
    cb = cos(-psia)
    sc = sin(-omegaa)
    cc = cos(-omegaa)
    sd = sin(chia)
    cd = cos(chia)

    # Compute elements of precession rotation matrix equivalent to
    # R3(chi_a) R1(-omega_a) R3(-psi_a) R1(epsilon_0).

    rot3 = array(((cd * cb - sb * sd * cc,
                   cd * sb * ca + sd * cc * cb * ca - sa * sd * sc,
                   cd * sb * sa + sd * cc * cb * sa + ca * sd * sc),
                  (-sd * cb - sb * cd * cc,
                   -sd * sb * ca + cd * cc * cb * ca - sa * cd * sc,
                   -sd * sb * sa + cd * cc * cb * sa + ca * cd * sc),
                  (sb * sc,
                   -sc * cb * ca - sa * cc,
                   -sc * cb * sa + cc * ca)))

    return rot3
import datetime as dt_module
import re
import sys
from collections import namedtuple
from datetime import date, datetime, timedelta
from numpy import (
    array, concatenate, cos, float_, int64, isnan, isinf, linspace,
    nan, ndarray, nonzero, pi, rollaxis, searchsorted, sin, where, zeros_like,

)
from time import strftime, struct_time
from ._compatibility import interp
from .constants import ASEC2RAD, B1950, DAY_S, T0, tau
from .curvelib import Splines, build_spline_given_ends
from .descriptorlib import reify
from .earthlib import sidereal_time, earth_rotation_angle
from .framelib import ICRS_to_J2000 as B
from .functions import (A, mxm, mxmxm, load_bundled_npy, rot_x, rot_y, rot_z,
                        _to_array, _reconcile)
from .nutationlib import (
    build_nutation_matrix, equation_of_the_equinoxes_complimentary_terms,
    iau2000a_radians, mean_obliquity,
)
from .precessionlib import compute_precession
        from astropy.time import Time

    def tt_strftime(self, format='%Y-%m-%d %H:%M:%S TT'):
        """Format TT with a datetime strftime() format string."""
        return self.ts._strftime(format, self.whole, self.tt_fraction)

    def tdb_strftime(self, format='%Y-%m-%d %H:%M:%S TDB'):
        """Format TDB with a datetime strftime() format string."""
        return self.ts._strftime(format, self.whole, self.tdb_fraction)

    def ut1_strftime(self, format='%Y-%m-%d %H:%M:%S UT1'):
        """Format UT1 with a datetime strftime() format string."""
        return self.ts._strftime(format, self.whole, self.ut1_fraction)

    # Convenient caching of several expensive functions of time.

    @reify
    def M(self):
        """3×3 rotation matrix: ICRS → equinox of this date."""

        # Compute N and P instead of asking for self.N and self.P to
        # avoid keeping copies of them, since Skyfield itself never uses
        # them again once M has been computed.  But we do check in case
        # a user has forced them to be built already.

        d = self.__dict__

        P = d.get('P')
        if P is None:
            P = self.precession_matrix()

        N = d.get('N')
        if N is None:
            N = self.nutation_matrix()

        return mxmxm(N, P, B)

    @reify
    def MT(self):
        """3×3 rotation matrix: equinox of this date → ICRS."""
        return rollaxis(self.M, 1)

    @reify
    def C(self):
        # Calculate the Equation of Origins in cycles
        eq_origins = (earth_rotation_angle(self.ut1) - self.gast / 24.0)
        R = rot_z(2 * pi * eq_origins)
        return mxm(R, self.M)

    @reify
    def CT(self):
        return rollaxis(self.C, 1)

    @reify
    def _nutation_angles_radians(self):
        # TODO: add psi and eps corrections support back in here, rather
        # than at points of use.
        return iau2000a_radians(self)

    def _nutation_angles(self, angles):
        # Sample code shared with early adopters suggested that setting
        # this attribute manually could avoid the expense of IAU 2000A,
        # so this setter continues to support the pattern.

        d_psi, d_eps = angles
        self._nutation_angles_radians = (
            d_psi / 1e7 * ASEC2RAD,
            d_eps / 1e7 * ASEC2RAD,
        )

    _nutation_angles = property(None, _nutation_angles)

    @reify
    def _mean_obliquity_radians(self):
        # Cached because it is used to compute both gast and N.
        return mean_obliquity(self.tdb) * ASEC2RAD

    # Conversion between timescales.

    @reify
    def J(self):
        """Return a floating point Julian year or array of years for this date.

        Julian years are convenient uniform periods of exactly 365.25
        days of Terrestrial Time, centered on 2000 January 1 12h TT =
        Julian year 2000.0.

        """
        return (self.whole - 1721045.0 + self.tt_fraction) / 365.25

    @reify
    def utc(self):
        """A tuple ``(year, month, day, hour, minute, seconds)`` in UTC."""
        utc = self._utc_tuple(0.0)
        return (array(utc).view(CalendarArray) if self.shape
                else CalendarTuple(*utc))

    @reify
    def tai_fraction(self):
        return self.tt_fraction - tt_minus_tai

    @reify
    def tdb_fraction(self):
        fr = self.tt_fraction
        return fr + tdb_minus_tt(self.whole, fr) / DAY_S

    @reify
    def ut1_fraction(self):
        return self.tt_fraction - self.delta_t / DAY_S

    @reify
    def delta_t(self):
        return self.ts.delta_t_function(self.tt)

    @reify
    def dut1(self):
        return 32.184 + self._leap_seconds() - self.delta_t

    @reify
    def gmst(self):
        """Greenwich Mean Sidereal Time (GMST) in hours."""

        return sidereal_time(self)