# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/nutationlib.py
# def iau2000a(jd_tt, fundamental_argument_terms=5, lunisolar_terms=687,
#              planetary_terms=687):
#     """Compute Earth nutation based on the IAU 2000A 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 value is either a float, or a NumPy array
#     with the same dimensions as the input argument.
# 
#     Supply smaller integer values for ``fundamental_argument_terms``,
#     ``lunisolar_terms``, and ``planetary_terms`` to trade off accuraccy
#     for speed.
# 
#     """
#     # Interval between fundamental epoch J2000.0 and given date.
# 
#     t = (jd_tt - T0) / 36525.0
# 
#     # Compute fundamental arguments from Simon et al. (1994), in radians.
# 
#     a = fundamental_arguments(t, fundamental_argument_terms)
# 
#     # ** Luni-solar nutation **
#     # Summation of luni-solar nutation series.
# 
#     cutoff = lunisolar_terms
#     arg = nals_t[:cutoff].dot(a).T
# 
#     sarg = sin(arg)
#     carg = cos(arg)
# 
#     dpsi = dot(sarg, lunisolar_longitude_coefficients[:cutoff,0])
#     dpsi += dot(sarg, lunisolar_longitude_coefficients[:cutoff,1]) * t
#     dpsi += dot(carg, lunisolar_longitude_coefficients[:cutoff,2])
# 
#     deps = dot(carg, lunisolar_obliquity_coefficients[:cutoff,0])
#     deps += dot(carg, lunisolar_obliquity_coefficients[:cutoff,1]) * t
#     deps += dot(sarg, lunisolar_obliquity_coefficients[:cutoff,2])
# 
#     # Compute and add in planetary components.
# 
#     if not planetary_terms:
#         return dpsi, deps
# 
#     if getattr(t, 'shape', ()) == ():
#         a = t * anomaly_coefficient + anomaly_constant
#     else:
#         a = (outer(anomaly_coefficient, t).T + anomaly_constant).T
#     a[-1] *= t
# 
#     cutoff = planetary_terms
#     arg = napl_t[:cutoff].dot(a).T
# 
#     sarg = sin(arg)
#     carg = cos(arg)
# 
#     dpsi += dot(sarg, nutation_coefficients_longitude[:cutoff,0])
#     dpsi += dot(carg, nutation_coefficients_longitude[:cutoff,1])
# 
#     deps += dot(sarg, nutation_coefficients_obliquity[:cutoff,0])
#     deps += dot(carg, nutation_coefficients_obliquity[:cutoff,1])
# 
#     return dpsi, deps
# 
# 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 matplotlib.pyplot as plt
from numpy import sqrt
from skyfield.api import load
from skyfield.nutationlib import iau2000a, iau2000b
"""Print a table of how accurate IAU 2000B is with fewer coefficients."""


ts = load.timescale()
#t = ts.utc(2000, 1, range(366 * 50))
t = ts.utc(1995, 1, range(366 * (2050 - 1995)))
print(t[0].utc_jpl())
print(t[-1].utc_jpl())


# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/nutationlib.py
# def iau2000a(jd_tt, fundamental_argument_terms=5, lunisolar_terms=687,
#              planetary_terms=687):
#     """Compute Earth nutation based on the IAU 2000A 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 value is either a float, or a NumPy array
#     with the same dimensions as the input argument.
# 
#     Supply smaller integer values for ``fundamental_argument_terms``,
#     ``lunisolar_terms``, and ``planetary_terms`` to trade off accuraccy
#     for speed.
# 
#     """
#     # Interval between fundamental epoch J2000.0 and given date.
# 
#     t = (jd_tt - T0) / 36525.0
# 
#     # Compute fundamental arguments from Simon et al. (1994), in radians.
# 
#     a = fundamental_arguments(t, fundamental_argument_terms)
# 
#     # ** Luni-solar nutation **
#     # Summation of luni-solar nutation series.
# 
#     cutoff = lunisolar_terms
#     arg = nals_t[:cutoff].dot(a).T
# 
#     sarg = sin(arg)
#     carg = cos(arg)
# 
#     dpsi = dot(sarg, lunisolar_longitude_coefficients[:cutoff,0])
#     dpsi += dot(sarg, lunisolar_longitude_coefficients[:cutoff,1]) * t
#     dpsi += dot(carg, lunisolar_longitude_coefficients[:cutoff,2])
# 
#     deps = dot(carg, lunisolar_obliquity_coefficients[:cutoff,0])
#     deps += dot(carg, lunisolar_obliquity_coefficients[:cutoff,1]) * t
#     deps += dot(sarg, lunisolar_obliquity_coefficients[:cutoff,2])
# 
#     # Compute and add in planetary components.
# 
#     if not planetary_terms:
#         return dpsi, deps
# 
#     if getattr(t, 'shape', ()) == ():
#         a = t * anomaly_coefficient + anomaly_constant
#     else:
#         a = (outer(anomaly_coefficient, t).T + anomaly_constant).T
#     a[-1] *= t
# 
#     cutoff = planetary_terms
#     arg = napl_t[:cutoff].dot(a).T
# 
#     sarg = sin(arg)
#     carg = cos(arg)
# 
#     dpsi += dot(sarg, nutation_coefficients_longitude[:cutoff,0])
#     dpsi += dot(carg, nutation_coefficients_longitude[:cutoff,1])
# 
#     deps += dot(sarg, nutation_coefficients_obliquity[:cutoff,0])
#     deps += dot(carg, nutation_coefficients_obliquity[:cutoff,1])
# 
#     return dpsi, deps
# 
# 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/nutation_table.py
import matplotlib.pyplot as plt
from numpy import sqrt
from skyfield.api import load
from skyfield.nutationlib import iau2000a, iau2000b

"""Print a table of how accurate IAU 2000B is with fewer coefficients."""


ts = load.timescale()
#t = ts.utc(2000, 1, range(366 * 50))
t = ts.utc(1995, 1, range(366 * (2050 - 1995)))
print(t[0].utc_jpl())
print(t[-1].utc_jpl())


dpsi_a, depsilon_a = iau2000a(t.tt)