# Path: skyfield/almanac.py
# def phase_angle(ephemeris, body, t):
# def fraction_illuminated(ephemeris, body, t):
# def seasons(ephemeris):
#     def season_at(t):
# def moon_phase(ephemeris, t):
# def moon_phases(ephemeris):
#     def moon_phase_at(t):
# def moon_nodes(ephemeris):
#     def moon_node_at(t):
# def oppositions_conjunctions(ephemeris, target):
#     def leading_or_trailing(t):
# def meridian_transits(ephemeris, target, topos):
#     def west_of_meridian_at(t):
# def sunrise_sunset(ephemeris, topos):
#     def is_sun_up_at(t):
# def dark_twilight_day(ephemeris, topos):
#     def is_it_dark_twilight_day_at(t):
# def risings_and_settings(ephemeris, target, topos,
#                          horizon_degrees=-34.0/60.0, radius_degrees=0): #?
#     def is_body_up_at(t):
# 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',
# ]
# MOON_PHASES = [
#     'New Moon',
#     'First Quarter',
#     'Full Moon',
#     'Last Quarter',
# ]
# MOON_NODES = [
#     'descending',
#     'ascending',
# ]
# CONJUNCTIONS = [
#     'conjunction',
#     'opposition',
# ]
# MERIDIAN_TRANSITS = ['Antimeridian transit', 'Meridian transit']
# TWILIGHTS = {
#     0: 'Night',
#     1: 'Astronomical twilight',
#     2: 'Nautical twilight',
#     3: 'Civil twilight',
#     4: 'Day',
# }
# 
# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/magnitudelib.py
# def planetary_magnitude(position):
#     """Given the position of a planet, return its visual magnitude.
# 
#     >>> from skyfield.api import load
#     >>> from skyfield.magnitudelib import planetary_magnitude
#     >>> ts = load.timescale()
#     >>> t = ts.utc(2020, 7, 31)
#     >>> eph = load('de421.bsp')
#     >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter'])
#     >>> print('%.2f' % planetary_magnitude(astrometric))
#     -2.73
# 
#     The formulae are from `Mallama and Hilton “Computing Apparent
#     Planetary Magnitude for the Astronomical Almanac” (2018)
#     <https://arxiv.org/pdf/1808.01973.pdf>`_.  Two of the formulae have
#     inherent limits:
# 
#     * Saturn’s magnitude is unknown and the function will return ``nan``
#       (the floating-point value “Not a Number”) if the “illumination
#       phase angle” — the angle of the vertex observer-Saturn-Sun —
#       exceeds 6.5°.
# 
#     * Neptune’s magnitude is unknown and will return ``nan`` if the
#       illumination phase angle exceeds 1.9° and the position's date is
#       before the year 2000.
# 
#     And one formula is not fully implemented (though contributions are
#     welcome!):
# 
#     * Skyfield does not compute which features on Mars are facing the
#       observer, which can introduce an error of ±0.06 magnitude.
# 
#     """
#     target = position.target
#     function = _FUNCTIONS.get(target)
#     if function is None:
#         name = _target_name(target)
#         raise ValueError('cannot compute the magnitude of target %s' % name)
# 
#     # Shamelessly treat the Sun as sitting at the Solar System Barycenter.
#     sun_to_observer = position.center_barycentric.xyz.au
#     observer_to_planet = position.xyz.au
#     sun_to_planet = sun_to_observer + observer_to_planet
# 
#     r = length_of(sun_to_planet)
#     delta = length_of(observer_to_planet)
#     ph_ang = angle_between(sun_to_planet, observer_to_planet) * RAD2DEG
# 
#     if function is _saturn_magnitude:
#         a = angle_between(_SATURN_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_SATURN_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _uranus_magnitude:
#         a = angle_between(_URANUS_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_URANUS_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _neptune_magnitude:
#         year = position.t.J
#         return function(r, delta, ph_ang, year)
# 
#     return function(r, delta, ph_ang)


import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from skyfield import almanac
from skyfield.api import load, wgs84
from skyfield.magnitudelib import planetary_magnitude

MONTH_NAMES = '0 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split()

# Figure out the times of sunset over our range of dates.

eph = load('de421.bsp')
earth, sun, venus = eph['earth'], eph['sun'], eph['venus']

# Path: skyfield/almanac.py
# def phase_angle(ephemeris, body, t):
# def fraction_illuminated(ephemeris, body, t):
# def seasons(ephemeris):
#     def season_at(t):
# def moon_phase(ephemeris, t):
# def moon_phases(ephemeris):
#     def moon_phase_at(t):
# def moon_nodes(ephemeris):
#     def moon_node_at(t):
# def oppositions_conjunctions(ephemeris, target):
#     def leading_or_trailing(t):
# def meridian_transits(ephemeris, target, topos):
#     def west_of_meridian_at(t):
# def sunrise_sunset(ephemeris, topos):
#     def is_sun_up_at(t):
# def dark_twilight_day(ephemeris, topos):
#     def is_it_dark_twilight_day_at(t):
# def risings_and_settings(ephemeris, target, topos,
#                          horizon_degrees=-34.0/60.0, radius_degrees=0): #?
#     def is_body_up_at(t):
# 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',
# ]
# MOON_PHASES = [
#     'New Moon',
#     'First Quarter',
#     'Full Moon',
#     'Last Quarter',
# ]
# MOON_NODES = [
#     'descending',
#     'ascending',
# ]
# CONJUNCTIONS = [
#     'conjunction',
#     'opposition',
# ]
# MERIDIAN_TRANSITS = ['Antimeridian transit', 'Meridian transit']
# TWILIGHTS = {
#     0: 'Night',
#     1: 'Astronomical twilight',
#     2: 'Nautical twilight',
#     3: 'Civil twilight',
#     4: 'Day',
# }
# 
# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/magnitudelib.py
# def planetary_magnitude(position):
#     """Given the position of a planet, return its visual magnitude.
# 
#     >>> from skyfield.api import load
#     >>> from skyfield.magnitudelib import planetary_magnitude
#     >>> ts = load.timescale()
#     >>> t = ts.utc(2020, 7, 31)
#     >>> eph = load('de421.bsp')
#     >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter'])
#     >>> print('%.2f' % planetary_magnitude(astrometric))
#     -2.73
# 
#     The formulae are from `Mallama and Hilton “Computing Apparent
#     Planetary Magnitude for the Astronomical Almanac” (2018)
#     <https://arxiv.org/pdf/1808.01973.pdf>`_.  Two of the formulae have
#     inherent limits:
# 
#     * Saturn’s magnitude is unknown and the function will return ``nan``
#       (the floating-point value “Not a Number”) if the “illumination
#       phase angle” — the angle of the vertex observer-Saturn-Sun —
#       exceeds 6.5°.
# 
#     * Neptune’s magnitude is unknown and will return ``nan`` if the
#       illumination phase angle exceeds 1.9° and the position's date is
#       before the year 2000.
# 
#     And one formula is not fully implemented (though contributions are
#     welcome!):
# 
#     * Skyfield does not compute which features on Mars are facing the
#       observer, which can introduce an error of ±0.06 magnitude.
# 
#     """
#     target = position.target
#     function = _FUNCTIONS.get(target)
#     if function is None:
#         name = _target_name(target)
#         raise ValueError('cannot compute the magnitude of target %s' % name)
# 
#     # Shamelessly treat the Sun as sitting at the Solar System Barycenter.
#     sun_to_observer = position.center_barycentric.xyz.au
#     observer_to_planet = position.xyz.au
#     sun_to_planet = sun_to_observer + observer_to_planet
# 
#     r = length_of(sun_to_planet)
#     delta = length_of(observer_to_planet)
#     ph_ang = angle_between(sun_to_planet, observer_to_planet) * RAD2DEG
# 
#     if function is _saturn_magnitude:
#         a = angle_between(_SATURN_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_SATURN_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _uranus_magnitude:
#         a = angle_between(_URANUS_POLE, sun_to_planet)
#         sun_sub_lat = a * RAD2DEG - 90.0
# 
#         a = angle_between(_URANUS_POLE, observer_to_planet)
#         observer_sub_lat = a * RAD2DEG - 90.0
# 
#         return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat)
# 
#     if function is _neptune_magnitude:
#         year = position.t.J
#         return function(r, delta, ph_ang, year)
# 
#     return function(r, delta, ph_ang)


# Path: examples/venus_evening_chart.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from skyfield import almanac
from skyfield.api import load, wgs84
from skyfield.magnitudelib import planetary_magnitude


MONTH_NAMES = '0 Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split()

# Figure out the times of sunset over our range of dates.

eph = load('de421.bsp')
earth, sun, venus = eph['earth'], eph['sun'], eph['venus']

observer = wgs84.latlon(+40.0, 0.0)