"""Formulae for specific earth behaviors and effects."""

from numpy import (abs, arcsin, arccos, arctan2, array, clip, cos, minimum,
                   nan_to_num, 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

earth_radius_au = ERAD / AU_M
one_minus_flattening = 1.0 - 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
one_minus_flattening_squared = one_minus_flattening * one_minus_flattening

def terra(latitude, longitude, elevation, gast):
    """Deprecated conversion from lat,lon,t -> GCRS; neglects polar motion."""
    zero = zeros_like(gast)
    sinphi = sin(latitude)
    cosphi = cos(latitude)
    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
    lat = arctan2(z, R)

    a = ERAD / AU_M
    f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
    e2 = 2.0*f - f*f
    i = 0
    C = 1.0
    while i < iterations:
        i += 1
        C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
        lat = arctan2(z + a * C * e2 * sin(lat), R)
    elevation_m = ((R / cos(lat)) - a * C) * AU_M
    return lat, lon, elevation_m

def compute_limb_angle(position_au, observer_au):
    """Determine the angle of an object above or below the Earth's limb.

    Given an object's GCRS `position_au` |xyz| vector and the position
    of an `observer_au` as a vector in the same coordinate system,
    return a tuple that provides `(limb_ang, nadir_ang)`:

    limb_angle
        Angle of observed object above (+) or below (-) limb in degrees.
    nadir_angle
        Nadir angle of observed object as a fraction of apparent radius
        of limb: <1.0 means below the limb, =1.0 means on the limb, and
        >1.0 means above the limb.

    """
    # Compute the distance to the object and the distance to the observer.

    disobj = sqrt(dots(position_au, position_au))
    disobs = sqrt(dots(observer_au, observer_au))

    # Compute apparent angular radius of Earth's limb.

    aprad = arcsin(minimum(earth_radius_au / disobs, 1.0))

    # Compute zenith distance of Earth's limb.

    zdlim = pi - aprad

    # Compute zenith distance of observed object.

    coszd = dots(position_au, observer_au) / (disobj * disobs)
    coszd = clip(coszd, -1.0, 1.0)
    zdobj = arccos(coszd)

    # Angle of object wrt limb is difference in zenith distances.

    limb_angle = (zdlim - zdobj) * RAD2DEG

    # Nadir angle of object as a fraction of angular radius of limb.

    nadir_angle = (pi - zdobj) / aprad

    return limb_angle, nadir_angle

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 refraction(alt_degrees, temperature_C, pressure_mbar):
    """Given an observed altitude, estimate atmospheric refraction, in degrees.

    Zero refraction is returned both for objects very near the zenith,
    as well as for objects more than one degree below the horizon.

    """
    r = 0.016667 / tan((alt_degrees + 7.31 / (alt_degrees + 4.4)) * DEG2RAD)
    d = r * (0.28 * pressure_mbar / (temperature_C + 273.0))
    return where((-1.0 <= alt_degrees) & (alt_degrees <= 89.9), d, 0.0)

def refract(alt_degrees, temperature_C, pressure_mbar):
    """Given an unrefracted `alt` determine where it will appear in the sky."""
    alt = alt_degrees
    while True:
        alt1 = alt
        alt = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
        converged = nan_to_num(abs(alt - alt1)).max() <= 3.0e-5
        if converged:
            break
    return alt