import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.time import Time
from six.moves.urllib.error import HTTPError

# Download and cache the IERS Bulletins A and B  using astropy's machinery
# (reminder: astroplan has its own function for this: `download_IERS_A`)
from astropy.utils.iers import (IERS_A, IERS_A_URL, IERS_B, IERS_B_URL,
                                FROM_IERS_B, IERS_Auto)
from astropy.utils.data import download_file
# Workaround until astropy/astropy#5194 is solved
iers_a = IERS_Auto.open()
iers_b = IERS_B.open(download_file(IERS_B_URL, cache=True))

# Get a range of times to plot from 1990-2022
time_range = Time("1990-01-01") + np.arange(0, 30, 0.2)*u.year

# Calculate the difference between UTC and UT1 at those times,
# allowing times "outside of the table"
DUT1_a, success_a = time_range.get_delta_ut1_utc(return_status=True,
                                                 iers_table=iers_a)
DUT1_b, success_b = time_range.get_delta_ut1_utc(return_status=True,
                                                 iers_table=iers_b)

# Compare input times to the times available in the table. For details, see
# https://github.com/astropy/astropy/blob/master/astropy/utils/iers/iers.py#L80
measurements_from_b = (success_b == FROM_IERS_B)

# Make a plot of the time difference
fig, ax = plt.subplots(figsize=(10,8))
ax.axhline(0, color='gray', ls='--', lw=2)

ax.plot_date(time_range.plot_date,
             DUT1_a, '-', lw=2, label='IERS Bulletin A + extrapolation')
ax.plot_date(time_range.plot_date[measurements_from_b],
             DUT1_b[measurements_from_b], 'r--', lw=2, label='IERS Bulletin B')
ax.set(xlabel='Year', ylabel='UT1-UTC [seconds]')
ax.legend(loc='upper right')
plt.show()