# Path: skyfield/constants.py
# 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
# 
# Path: skyfield/functions.py
# 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)


from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos,
                   minimum, pi, sin, sqrt, tan, where, zeros_like)
from .constants import (AU_M, ANGVEL, DAY_S, DEG2RAD, ERAD,
                        IERS_2010_INVERSE_EARTH_FLATTENING, RAD2DEG, T0, tau)
from .functions import dots
    c = 1.0 / sqrt(cosphi * cosphi +
                   sinphi * sinphi * one_minus_flattening_squared)
    s = one_minus_flattening_squared * c
    ach = earth_radius_au * c + elevation
    ash = earth_radius_au * s + elevation

    # Compute local sidereal time factors at the observer's longitude.

    stlocl = 15.0 * DEG2RAD * gast + longitude
    sinst = sin(stlocl)
    cosst = cos(stlocl)

    # Compute position vector components in au.

    ac = ach * cosphi
    acsst = ac * sinst
    accst = ac * cosst
    pos = array((accst, acsst, zero + ash * sinphi))

    # Compute velocity vector components in au/day.

    vel = ANGVEL * DAY_S * array((-acsst, accst, zero))

    return pos, vel

def reverse_terra(xyz_au, gast, iterations=3):
    """Deprecated conversion from GCRS -> lat,lon,t; neglects polar motion."""
    x, y, z = xyz_au
    R = sqrt(x*x + y*y)


# Path: skyfield/constants.py
# 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
# 
# Path: skyfield/functions.py
# 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)


# Path: skyfield/earthlib.py
from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos,
                   minimum, pi, sin, sqrt, tan, where, zeros_like)
from .constants import (AU_M, ANGVEL, DAY_S, DEG2RAD, ERAD,
                        IERS_2010_INVERSE_EARTH_FLATTENING, RAD2DEG, T0, tau)
from .functions import dots

    c = 1.0 / sqrt(cosphi * cosphi +
                   sinphi * sinphi * one_minus_flattening_squared)
    s = one_minus_flattening_squared * c
    ach = earth_radius_au * c + elevation
    ash = earth_radius_au * s + elevation

    # Compute local sidereal time factors at the observer's longitude.

    stlocl = 15.0 * DEG2RAD * gast + longitude
    sinst = sin(stlocl)
    cosst = cos(stlocl)

    # Compute position vector components in au.

    ac = ach * cosphi
    acsst = ac * sinst
    accst = ac * cosst
    pos = array((accst, acsst, zero + ash * sinphi))

    # Compute velocity vector components in au/day.

    vel = ANGVEL * DAY_S * array((-acsst, accst, zero))

    return pos, vel

def reverse_terra(xyz_au, gast, iterations=3):
    """Deprecated conversion from GCRS -> lat,lon,t; neglects polar motion."""
    x, y, z = xyz_au
    R = sqrt(x*x + y*y)


    lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi