# -*- coding: utf-8 -*-
"""Routines to solve for circumstances like sunrise, sunset, and moon phase."""

from __future__ import print_function, division

from numpy import cos, zeros_like
from .constants import pi, tau
from .framelib import ecliptic_frame
from .searchlib import find_discrete
from .nutationlib import iau2000b_radians
from .units import Angle

# Not only to support historic code but also for future convenience, let
# folks import the search routine alongside the almanac routines.
find_discrete

# Simple facts.

def phase_angle(ephemeris, body, t):
    """
    .. deprecated:: 1.42
       Use the :meth:`~skyfield.positionlib.ICRF.phase_angle()` position
       method instead.
    """
    p = ephemeris['earth'].at(t).observe(ephemeris[body])
    return p.phase_angle(ephemeris['sun'])

def fraction_illuminated(ephemeris, body, t):
    """
    .. deprecated:: 1.42
       Use the :meth:`~skyfield.positionlib.ICRF.fraction_illuminated()`
       position method instead.
    """
    a = phase_angle(ephemeris, body, t).radians
    return 0.5 * (1.0 + cos(a))

# Discrete circumstances to search.

SEASONS = [
    'Spring',
    'Summer',
    'Autumn',
    'Winter',
]

SEASON_EVENTS = [
    'Vernal Equinox',
    'Summer Solstice',
    'Autumnal Equinox',
    'Winter Solstice',
]

SEASON_EVENTS_NEUTRAL = [
    'March Equinox',
    'June Solstice',
    'September Equinox',
    'December Solstice',
]

def seasons(ephemeris):
    """Build a function of time that returns the quarter of the year.

    The function that this returns will expect a single argument that is
    a :class:`~skyfield.timelib.Time` and will return 0 through 3 for
    the seasons Spring, Summer, Autumn, and Winter.

    """
    earth = ephemeris['earth']
    sun = ephemeris['sun']

    def season_at(t):
        """Return season 0 (Spring) through 3 (Winter) at time `t`."""
        t._nutation_angles_radians = iau2000b_radians(t)
        e = earth.at(t)
        _, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
        return (slon.radians // (tau / 4) % 4).astype(int)

    season_at.step_days = 90.0
    return season_at

MOON_PHASES = [
    'New Moon',
    'First Quarter',
    'Full Moon',
    'Last Quarter',
]

def moon_phase(ephemeris, t):
    """Return the Moon phase 0°–360° at time ``t``, where 180° is Full Moon.

    More precisely: this returns an :class:`~skyfield.units.Angle`
    giving the difference between the geocentric apparent ecliptic
    longitudes of the Moon and Sun, constrained to the interval 0°–360°
    (0–𝜏 radians) where 0° is New Moon and 180° is Full Moon.

    """
    e = ephemeris['earth'].at(t)
    moon, sun = ephemeris['moon'], ephemeris['sun']
    _, mlon, _ = e.observe(moon).apparent().frame_latlon(ecliptic_frame)
    _, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
    return Angle(radians=(mlon.radians - slon.radians) % tau)

def moon_phases(ephemeris):
    """Build a function of time that returns the moon phase 0 through 3.

    The function that this returns will expect a single argument that is
    a :class:`~skyfield.timelib.Time` and will return the phase of the
    moon as an integer.  See the accompanying array ``MOON_PHASES`` if
    you want to give string names to each phase.

    """
    earth = ephemeris['earth']
    moon = ephemeris['moon']
    sun = ephemeris['sun']

    def moon_phase_at(t):
        """Return the phase of the moon 0 through 3 at time `t`."""
        t._nutation_angles_radians = iau2000b_radians(t)
        e = earth.at(t)
        _, mlon, _ = e.observe(moon).apparent().frame_latlon(ecliptic_frame)
        _, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
        return ((mlon.radians - slon.radians) // (tau / 4) % 4).astype(int)

    moon_phase_at.step_days = 7.0  # one lunar phase per week
    return moon_phase_at

MOON_NODES = [
    'descending',
    'ascending',
]

def moon_nodes(ephemeris):
    """Build a function of time that identifies lunar nodes.

    This returns a function taking a :class:`~skyfield.timelib.Time` and
    returning ``True`` if the Moon is above the ecliptic else ``False``.
    See :ref:`lunar-nodes` for how to use this routine.

    """
    earth = ephemeris['earth']
    moon = ephemeris['moon']

    def moon_node_at(t):
        """Return the phase of the moon 0 through 3 at time `t`."""
        e = earth.at(t)
        lat, _, _ = e.observe(moon).apparent().frame_latlon(ecliptic_frame)
        return lat.radians > 0.0

    moon_node_at.step_days = 12.0  # 2000-2050: closest nodes 12.38 days apart
    return moon_node_at

CONJUNCTIONS = [
    'conjunction',
    'opposition',
]

def oppositions_conjunctions(ephemeris, target):
    """Build a function to find oppositions and conjunctions with the Sun.

    See :ref:`oppositions-conjunctions` for how to call this routine and
    interpret the results.

    """
    earth_at = ephemeris['earth'].at
    sun = ephemeris['sun']

    def leading_or_trailing(t):
        """Return whether the target is east or west of the Sun."""
        e = earth_at(t)
        _, slon, _ = e.observe(sun).apparent().frame_latlon(ecliptic_frame)
        _, tlon, _ = e.observe(target).apparent().frame_latlon(ecliptic_frame)
        return ((slon.radians - tlon.radians) / pi % 2.0).astype('int8')

    if target.target == 301:
        leading_or_trailing.step_days = 14  # Moon
    else:
        leading_or_trailing.step_days = 40  # Mercury (the fastest planet)
    return leading_or_trailing

MERIDIAN_TRANSITS = ['Antimeridian transit', 'Meridian transit']

def meridian_transits(ephemeris, target, topos):
    """Build a function of time for finding when a body transits the meridian.

    The returned function accepts a :class:`~skyfield.timelib.Time`
    argument and returns ``True`` if the ``target`` body is west of the
    observer’s meridian at that time, and otherwise returns ``False.``
    See :ref:`transits` for how to use this to search for a body’s
    meridian transits and antimeridian transits.

    """
    topos_at = (ephemeris['earth'] + topos).at

    def west_of_meridian_at(t):
        """Return `True` if the target is west of the observer’s meridian."""
        t._nutation_angles_radians = iau2000b_radians(t)
        # TODO: should we work to avoid computing Topos position twice?
        # We could grab its hidden GCRS vector and do the trig ourselves.
        # Or there might be something clever we can do with the two raw
        # vectors, skipping the cost of computing spherical coordinates.
        ra1, _, _ = topos.at(t).radec(epoch='date')
        ra2, _, _ = topos_at(t).observe(target).apparent().radec(epoch='date')
        return (ra1.radians - ra2.radians) % tau < pi

    west_of_meridian_at.step_days = 0.4  # twice a day
    return west_of_meridian_at

def sunrise_sunset(ephemeris, topos):
    """Build a function of time that returns whether the Sun is up.

    The function that is returned will expect a single argument that is
    a :class:`~skyfield.timelib.Time`, and will return ``True`` if the
    sun is up, else ``False``.

    Skyfield uses the same definition as the United States Naval
    Observatory: the Sun is up when its center is 0.8333 degrees below
    the horizon, which accounts for both its apparent radius of around
    16 arcminutes and also for the 34 arcminutes by which atmospheric
    refraction on average lifts the image of the Sun.

    If you need to provide a custom value for refraction, adjust the
    estimate of the Sun’s radius, or account for a vantage point above
    the Earth’s surface, see :ref:`risings-and-settings` to learn about
    the more versatile :func:`~skyfield.almanac.risings_and_settings()`
    routine.

    """
    sun = ephemeris['sun']
    topos_at = (ephemeris['earth'] + topos).at

    def is_sun_up_at(t):
        """Return `True` if the sun has risen by time `t`.

        The Sun has risen if its altitude above the horizon is greater
        than -0.8333 degrees.

        """
        t._nutation_angles_radians = iau2000b_radians(t)
        return topos_at(t).observe(sun).apparent().altaz()[0].degrees >= -0.8333

    is_sun_up_at.step_days = 0.04  # catch days at least an hour long
    return is_sun_up_at

TWILIGHTS = {
    0: 'Night',
    1: 'Astronomical twilight',
    2: 'Nautical twilight',
    3: 'Civil twilight',
    4: 'Day',
}

def dark_twilight_day(ephemeris, topos):
    """Build a function of time returning whether it is dark, twilight, or day.

    The function that this returns will expect a single argument that is
    a :class:`~skyfield.timelib.Time` and will return:

    | 0 — Dark of night.
    | 1 — Astronomical twilight.
    | 2 — Nautical twilight.
    | 3 — Civil twilight.
    | 4 — Sun is up.

    """
    sun = ephemeris['sun']
    topos_at = (ephemeris['earth'] + topos).at

    def is_it_dark_twilight_day_at(t):
        """Return whether the Sun is up, down, or whether there is twilight."""
        t._nutation_angles_radians = iau2000b_radians(t)
        degrees = topos_at(t).observe(sun).apparent().altaz()[0].degrees
        r = zeros_like(degrees, int)
        r[degrees >= -18.0] = 1
        r[degrees >= -12.0] = 2
        r[degrees >= -6.0] = 3
        r[degrees >= -0.8333] = 4
        return r

    is_it_dark_twilight_day_at.step_days = 0.04  # catch days at least an hour long
    return is_it_dark_twilight_day_at

def risings_and_settings(ephemeris, target, topos,
                         horizon_degrees=-34.0/60.0, radius_degrees=0): #?
    """Build a function of time that returns whether a body is up.

    This returns a function taking a :class:`~skyfield.timelib.Time`
    argument returning ``True`` if the body’s altazimuth altitude angle
    plus ``radius_degrees`` is greater than ``horizon_degrees``, else
    ``False``.  See :ref:`risings-and-settings` to learn about how to
    search for risings and settings, and to see more about using the
    parameters ``horizon_degrees`` and ``radius_degrees``.

    """
    topos_at = (ephemeris['earth'] + topos).at
    h = horizon_degrees - radius_degrees

    def is_body_up_at(t):
        """Return `True` if the target has risen by time `t`."""
        t._nutation_angles_radians = iau2000b_radians(t)
        return topos_at(t).observe(target).apparent().altaz()[0].degrees > h

    is_body_up_at.step_days = 0.25
    return is_body_up_at