# 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/nutationlib.py
# def iau2000b(jd_tt):
#     """Compute Earth nutation based on the faster IAU 2000B nutation model.
# 
#     `jd_tt` - Terrestrial Time: Julian date float, or NumPy array of floats
# 
#     Returns a tuple ``(delta_psi, delta_epsilon)`` measured in tenths of
#     a micro-arcsecond.  Each is either a float, or a NumPy array with
#     the same dimensions as the input argument.  The result will not take
#     as long to compute as the full IAU 2000A series, but should still
#     agree with ``iau2000a()`` to within a milliarcsecond between the
#     years 1995 and 2020.
# 
#     """
#     dpsi, deps = iau2000a(jd_tt, 2, 77, 0)
#     dpsi += -0.000135e7
#     deps +=  0.000388e7
#     return dpsi, deps


import numpy as np
from skyfield import almanac, api
from skyfield.nutationlib import iau2000b
from time import time
"""

ts = api.load.timescale()
eph = api.load('de421.bsp')
sat = api.EarthSatellite(
    '1 25544U 98067A   20150.54009403  .00000886  00000-0  23936-4 0  9992',
    '2 25544  51.6447  80.7632 0002782   3.7438 167.9307 15.49399736229124',
)

topos = api.Topos('18.3368 N', '64.7281 W')
t0 = ts.utc(2020, 6, 5)
t1 = ts.utc(2020, 6, 6)

half_second = 0.5 / 24.0 / 3600.0


T0 = time()
t, y = sat.find_events(topos, t0, t1, altitude_degrees=10.0)
print(time() - T0, 'seconds to find events:')

for ti, yi in zip(t, y):
    print(' ', ti.utc_jpl(), yi)

offset = 0
assert list(y[offset:offset+3]) == [0, 1, 2]

pass_start, pass_end = t[offset+0], t[offset+2]

def f(t):
    length = t.tt.shape[0]

# 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/nutationlib.py
# def iau2000b(jd_tt):
#     """Compute Earth nutation based on the faster IAU 2000B nutation model.
# 
#     `jd_tt` - Terrestrial Time: Julian date float, or NumPy array of floats
# 
#     Returns a tuple ``(delta_psi, delta_epsilon)`` measured in tenths of
#     a micro-arcsecond.  Each is either a float, or a NumPy array with
#     the same dimensions as the input argument.  The result will not take
#     as long to compute as the full IAU 2000A series, but should still
#     agree with ``iau2000a()`` to within a milliarcsecond between the
#     years 1995 and 2020.
# 
#     """
#     dpsi, deps = iau2000a(jd_tt, 2, 77, 0)
#     dpsi += -0.000135e7
#     deps +=  0.000388e7
#     return dpsi, deps


# Path: design/satellite_passes.py
import numpy as np
from skyfield import almanac, api
from skyfield.nutationlib import iau2000b
from time import time

"""

ts = api.load.timescale()
eph = api.load('de421.bsp')
sat = api.EarthSatellite(
    '1 25544U 98067A   20150.54009403  .00000886  00000-0  23936-4 0  9992',
    '2 25544  51.6447  80.7632 0002782   3.7438 167.9307 15.49399736229124',
)

topos = api.Topos('18.3368 N', '64.7281 W')
t0 = ts.utc(2020, 6, 5)
t1 = ts.utc(2020, 6, 6)

half_second = 0.5 / 24.0 / 3600.0


T0 = time()
t, y = sat.find_events(topos, t0, t1, altitude_degrees=10.0)
print(time() - T0, 'seconds to find events:')

for ti, yi in zip(t, y):
    print(' ', ti.utc_jpl(), yi)

offset = 0
assert list(y[offset:offset+3]) == [0, 1, 2]

pass_start, pass_end = t[offset+0], t[offset+2]

def f(t):
    length = t.tt.shape[0]

    t._nutation_angles = iau2000b(t.tt)