For #452, switch Skyfield to IERS data for ∆T

builders/build_arrays.py

@@ -26,9 +26,9 @@ def main(argv):
     ]})
 
     f = load.open(iers.FINALS_URL)
-    mjd, dut1 = iers.parse_dut1_from_finals_all(f)
+    mjd_utc, dut1 = iers.parse_dut1_from_finals_all(f)
     delta_t_recent, leap_dates, leap_offsets = (
-        iers.build_timescale_arrays(mjd, dut1)
+        iers.build_timescale_arrays(mjd_utc, dut1)
     )
     savez_compressed(
         'skyfield/data/iers',

design/delta_t.py

@@ -36,13 +36,13 @@ def main(argv):
     #draw_plot_comparing_USNO_and_IERS_data()
 
     f = load.open('finals.all')
-    mjd, dut1 = iers.parse_dut1_from_finals_all(f)
+    mjd_utc, dut1 = iers.parse_dut1_from_finals_all(f)
     delta_t_recent, leap_dates, leap_offsets = (
-        iers.build_timescale_arrays(mjd, dut1)
+        iers.build_timescale_arrays(mjd_utc, dut1)
     )
 
     ts = load.timescale()
-    jd = mjd + 2400000.5
+    jd = mjd_utc + 2400000.5
     t = ts.tt_jd(jd)
 
     fig, ax = plt.subplots()

skyfield/data/iers.py

@@ -2,31 +2,36 @@
 
 See:
 https://datacenter.iers.org/eop.php
-ftp://ftp.iers.org/products/eop/rapid/standard/
+ftp://cddis.gsfc.nasa.gov/pub/products/iers/readme.finals2000A
 
 """
 import numpy as np
 import re
 
+from ..constants import DAY_S
+
 FINALS_URL = 'ftp://ftp.iers.org/products/eop/rapid/standard/finals2000A.all'
 _DUT1 = re.compile(b'^......(.........) ' + b'.' * 42 + b'(.\d........)', re.M)
 inf = float('inf')
 
 def parse_dut1_from_finals_all(f):
     data = np.fromregex(f, _DUT1, [
-        ('mjd', np.float32),
+        ('mjd_utc', np.float32),
         ('dut1', np.float32),
     ])
-    return data['mjd'], data['dut1']
+    return data['mjd_utc'], data['dut1']
 
-def build_timescale_arrays(mjd, dut1):
+def build_timescale_arrays(mjd_utc, dut1):
     big_jumps = np.diff(dut1) > 0.9
     leap_second_mask = np.concatenate([[False], big_jumps])
-    delta_t = np.cumsum(leap_second_mask) - dut1 + 32.184 + 12.0
-    delta_t_recent = np.array([mjd + 2400000.5, delta_t])
+    tt_minus_utc = np.cumsum(leap_second_mask) + 32.184 + 12.0
+
+    tt_jd = mjd_utc + tt_minus_utc / DAY_S + 2400000.5
+    delta_t = tt_minus_utc - dut1
+    delta_t_recent = np.array([tt_jd, delta_t])
 
     leap_dates = 2400000.5 + np.concatenate([
-        [-inf], [41317.0, 41499.0, 41683.0], mjd[leap_second_mask], [inf],
+        [-inf], [41317.0, 41499.0, 41683.0], mjd_utc[leap_second_mask], [inf],
     ])
     leap_offsets = np.arange(8.0, len(leap_dates) + 8.0)
     leap_offsets[0] = leap_offsets[1] = 10.0

skyfield/documentation/almanac.rst

@@ -391,7 +391,7 @@ The difference in rising time can be a minute or more:
 
 .. testoutput::
 
-    2020-02-01 12:46:28Z Limb of the Sun crests the horizon
+    2020-02-01 12:46:27Z Limb of the Sun crests the horizon
     2020-02-01 12:47:53Z Center of the Sun reaches the horizon
 
 Elevated vantage points

skyfield/documentation/earth-satellites.rst

@@ -365,7 +365,7 @@ just as you did for the position measured from the Earth’s center:
 
 .. testoutput::
 
-    [ 331.61883722  392.18468536 1049.76011302]
+    [ 331.61874618  392.18485533 1049.76011316]
 
 But the most popular approach is to ask the topocentric position
 for its altitude and azimuth coordinates,

skyfield/documentation/examples.rst

@@ -57,7 +57,7 @@ or how early I need to rise to see the morning sky:
 
     2020-04-19 05:09   Astronomical twilight starts
     2020-04-19 05:46   Nautical twilight starts
-    2020-04-19 06:21   Civil twilight starts
+    2020-04-19 06:20   Civil twilight starts
     2020-04-19 06:49   Day starts
     2020-04-19 20:20   Civil twilight starts
     2020-04-19 20:48   Nautical twilight starts

skyfield/documentation/files.rst

@@ -3,61 +3,38 @@
  Downloading and Using Data Files
 ==================================
 
-Your Skyfield programs will typically download two kinds of data file.
-
-First, Skyfield will need up-to-date tables about time —
-files providing recently measured values for ∆T,
-future predictions of ∆T, and a table of recent leap seconds.
-See :ref:`downloading-timescale-files` to learn more about these files.
-
-Second, Skyfield will need data
-about the objects that you want to observe:
-planets, stars, or Earth satellites.
-The rest of the documentation
-will introduce you to each of these kinds of data source.
-
-The result will be a program that, when run for the first time,
-downloads several data files before it can perform useful operations.
-If the program’s output is a terminal window,
-then it will display progress bars as each file is downloaded:
-
-.. testsetup::
-
-   import os
-   from skyfield.api import Loader
-
-   load = Loader('.')
-   ts = load.timescale(builtin=False)
-   planets = load('de421.bsp')
+The first time you run a Skyfield program,
+it will typically download one or more data files from the Internet
+that help it compute the positions of planets or satellites —
+one file for each call the program makes to Skyfield’s ``load()`` routine.
+If the program is attached to a terminal,
+then a simple progress bar will be displayed
+as Skyfield downloads each file.
 
 ::
 
    from skyfield.api import load
-
-   ts = load.timescale(builtin=False)
    planets = load('de421.bsp')
    print('Ready')
 
 ::
 
-   [#################################] 100% deltat.data
-   [#################################] 100% deltat.preds
-   [#################################] 100% Leap_Second.dat
    [#################################] 100% de421.bsp
    Ready
 
 The second time you run the program, however,
-it will find the data files already sitting in the current directory
-and can start up without needing to access the network:
+the program will find the data files
+already sitting in the current directory.
+In that case, the program will run without needing access to the Internet:
 
 ::
 
    Ready
 
-Most programs will run just fine using the default ``load()`` object
+Most programs will run just fine using the default ``load()`` routine
 provided in the ``skyfield.api`` module.
 But other programs may want to build their own loader
-so that they have the chance to specify non-default behavior.
+so they have the chance to override its default behaviors.
 
 Specifying the download directory
 =================================
@@ -96,4 +73,5 @@ by building a `Loader` whose verbosity is set to false.
 
 .. testcode::
 
+   from skyfield.api import Loader
    load = Loader('~/skyfield-data', verbose=False)

skyfield/documentation/index.rst

@@ -71,7 +71,7 @@ on the Earth’s surface:
 .. testoutput::
 
     25deg 27' 54.0"
-    101deg 33' 44.0"
+    101deg 33' 44.1"
 
 Skyfield does not depend on the `AstroPy`_ project
 or its compiled libraries.

skyfield/iokit.py

@@ -11,12 +11,14 @@
 
 import certifi
 
+from .data import iers
+from .functions import load_bundled_npy
 from .io_timescale import (
-     _build_timescale, parse_deltat_data,
-    parse_deltat_preds, parse_leap_seconds,
+    parse_deltat_data, parse_deltat_preds, parse_leap_seconds,
 )
 from .jpllib import SpiceKernel
 from .sgp4lib import EarthSatellite
+from .timelib import Timescale
 
 try:
     from io import BytesIO
@@ -279,24 +281,31 @@ def open(self, url, mode='rb', reload=False, filename=None, backup=False):
         return open(path, mode)
 
     def timescale(self, delta_t=None, builtin=True):
-        """Open or download three time scale files, returning a `Timescale`.
+        """Return a `Timescale` built using official Earth rotation data.
 
         ``delta_t`` — Lets you override the standard ∆T tables by
         providing your own ∆T offset in seconds.  For details, see
         :ref:`custom-delta-t`.
 
-        ``builtin`` — Set this to ``False`` to download new copies of
-        the three official files that supply recent ∆T measurements and
-        the current schedule of UTC leap seconds.  By default, Skyfield
-        uses copies of the three files that it ships with and installs
-        alongside its code.  For details, see
+        ``builtin`` — By default, Skyfield uses ∆T and leap second
+        tables that it carries internally.  Set this option to ``False``
+        to download up-to-date data (about 3.3 MB) directly from the
+        International Earth Rotation Service instead.  For details, see
         :ref:`downloading-timescale-files`.
 
         """
-        deltat_data = self('deltat.data', builtin=builtin)
-        deltat_preds = self('deltat.preds', builtin=builtin)
-        _, leap_second_dat = self('Leap_Second.dat', builtin=builtin)
-        ts = _build_timescale(deltat_data, deltat_preds, leap_second_dat)
+        if builtin:
+            arrays = load_bundled_npy('iers.npz')
+            delta_t_recent = arrays['delta_t_recent']
+            leap_dates = arrays['leap_dates']
+            leap_offsets = arrays['leap_offsets']
+        else:
+            with self.open(iers.FINALS_URL) as f:
+                mjd_utc, dut1 = iers.parse_dut1_from_finals_all(f)
+            delta_t_recent, leap_dates, leap_offsets = (
+                iers.build_timescale_arrays(mjd_utc, dut1))
+
+        ts = Timescale(delta_t_recent, leap_dates, leap_offsets)
 
         if delta_t is not None:
             ts.delta_t_table = np.array(((-1e99, 1e99), (delta_t, delta_t)))

skyfield/tests/test_positions.py

@@ -5,7 +5,7 @@
 from skyfield.functions import length_of, mxv, rot_z
 from skyfield.positionlib import ICRF, ITRF_to_GCRS2, _GIGAPARSEC_AU
 from skyfield.starlib import Star
-from .fixes import IS_32_BIT, low_precision_ERA
+from .fixes import low_precision_ERA
 
 def test_subtraction():
     p0 = ICRF((10,20,30), (40,50,60))
@@ -157,7 +157,7 @@ def test_velocity_in_ITRF_to_GCRS2():
     relative_error = (length_of(actual_motion - predicted_motion)
                       / length_of(actual_motion))
 
-    acceptable_error = 4e-12 if IS_32_BIT else 2e-12
+    acceptable_error = 4e-12
     assert relative_error < acceptable_error
 
 # Test that the CIRS coordinate of the TIO is consistent with the Earth Rotation Angle

skyfield/tests/test_timelib.py

@@ -137,14 +137,14 @@ def test_strftime_with_microseconds(ts):
 
     t = ts.tt(2020, 9, 12)
     assert t.utc_strftime('%Y %S %f') == '2020 50 816000'
-    assert t.ut1_strftime('%Y %S %f') == '2020 49 776521'
+    assert t.ut1_strftime('%Y %S %f') == '2020 50 638465'
     assert t.tai_strftime('%Y %S %f') == '2020 27 816000'
     assert t.tt_strftime('%Y %S %f') == '2020 00 000000'
     assert t.tdb_strftime('%Y %S %f') == '2020 59 998446'
 
     t = ts.tt(2020, 9, [12, 12])
     assert t.utc_strftime('%Y %S %f') == ['2020 50 816000'] * 2
-    assert t.ut1_strftime('%Y %S %f') == ['2020 49 776521'] * 2
+    assert t.ut1_strftime('%Y %S %f') == ['2020 50 638465'] * 2
     assert t.tai_strftime('%Y %S %f') == ['2020 27 816000'] * 2
     assert t.tt_strftime('%Y %S %f') == ['2020 00 000000'] * 2
     assert t.tdb_strftime('%Y %S %f') == ['2020 59 998446'] * 2
@@ -436,9 +436,10 @@ def test_leap_second(ts):
     assert ts.tai(jd=t5).utc_iso() == '1974-01-01T00:00:01Z'
 
 def test_delta_t(ts):
-    # Check delta_t calculation around year 2000/1/1 (from IERS tables this is 63.8285)
+    # The IERS "finals2000A.all" for 2000 Jan 1 gives DUT1 = 0.3554779,
+    # and 0.3554779 - 0.184 - 1.0 = -0.8285221.
     t = ts.utc(2000, 1, 1, 0, 0, 0)
-    assert abs(t.delta_t - 63.8285) < 1e-5
+    assert abs(t.delta_t - 63.8285221) < 1e-9
 
     # Check historic value. Compare to the table in Morrison and
     # Stephenson 2004, the tolerance is 2 sigma