From 787d6c1b54851bd9a82025c7b9a270a94d112bc4 Mon Sep 17 00:00:00 2001
From: Brandon Rhodes <brandon@rhodesmill.org>
Date: Sun, 20 May 2018 14:13:06 -0400
Subject: [PATCH] Allow time objects to be constructed from UT1

Fixes #91.
---
 CHANGELOG.rst                  |  3 +++
 skyfield/tests/test_timelib.py | 17 ++++++++++++++-
 skyfield/timelib.py            | 39 ++++++++++++++++++++++++++++++++++
 3 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3bade737d..e446c4260 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -19,6 +19,9 @@ Changelog
 * Star coordinates can now be offered for any epoch, not just J2000.
   `#166 <https://github.com/skyfielders/python-skyfield/issues/166>`_
 
+* You can now create a time object given the UT1 date.
+  `#91 <https://github.com/skyfielders/python-skyfield/issues/91>`_
+
 1.3 — 2018 April 15
 -------------------
 
diff --git a/skyfield/tests/test_timelib.py b/skyfield/tests/test_timelib.py
index 86b9778fe..9c711c53e 100644
--- a/skyfield/tests/test_timelib.py
+++ b/skyfield/tests/test_timelib.py
@@ -9,7 +9,7 @@
 one_second = 1.0 / DAY_S
 epsilon = one_second * 42.0e-6  # 20.1e-6 is theoretical best precision
 
-time_parameter = ['tai', 'tt', 'tdb']
+time_parameter = ['tai', 'tt', 'tdb', 'ut1']
 time_value = [(1973, 1, 18, 1, 35, 37.5), 2441700.56640625]
 
 def ts():
@@ -66,6 +66,21 @@ def test_building_time_from_list_of_utc_datetimes(ts):
         2442046.5, 2442047.5, 2442048.5, 2442049.5, 2442050.5, 2442051.5,
         ]).all()
 
+def test_converting_ut1_to_tt(ts):
+    ten_thousand_years = 365 * 10000
+
+    jd = api.T0 - ten_thousand_years
+    t = ts.ut1(jd=jd)
+    del t.ut1                   # force re-computation of UT1
+    print(jd - t.ut1)
+    assert abs(jd - t.ut1) < 1e-10
+
+    jd = api.T0 + ten_thousand_years
+    t = ts.ut1(jd=jd)
+    del t.ut1                   # force re-computation of UT1
+    print(jd - t.ut1)
+    assert abs(jd - t.ut1) < 1e-10
+
 def test_indexing_time(ts):
     t = ts.utc(1974, 10, range(1, 6))
     assert t.shape == (5,)
diff --git a/skyfield/timelib.py b/skyfield/timelib.py
index 6ff6b542b..3da9a17ae 100644
--- a/skyfield/timelib.py
+++ b/skyfield/timelib.py
@@ -197,6 +197,45 @@ def tdb(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0,
         t.tdb = tdb
         return t
 
+    def ut1(self, year=None, month=1, day=1, hour=0, minute=0, second=0.0,
+            jd=None):
+        """Return the Time corresponding to a specific moment in UT1.
+
+        You can supply the Universal Time (UT1) by providing either a
+        proleptic Gregorian calendar date or a raw Julian Date float.
+        The following two method calls are equivalent::
+
+            timescale.ut1(2014, 1, 18, 1, 35, 37.5)
+            timescale.ut1(jd=2456675.56640625)
+
+        """
+        if jd is not None:
+            ut1 = jd
+        else:
+            ut1 = julian_date(
+                _to_array(year), _to_array(month), _to_array(day),
+                _to_array(hour), _to_array(minute), _to_array(second),
+            )
+        ut1 = _to_array(ut1)
+
+        # Estimate TT = UT1, to get a rough Delta T estimate.
+        tt_approx = ut1
+        delta_t_approx = interpolate_delta_t(self.delta_t_table, tt_approx)
+
+        # Use the rough Delta T to make a much better estimate of TT,
+        # then generate an even better Delta T.
+        tt_approx = ut1 + delta_t_approx / DAY_S
+        delta_t_approx = interpolate_delta_t(self.delta_t_table, tt_approx)
+
+        # We can now estimate TT with an error of < 1e-9 seconds within
+        # 10 centuries of either side of the present; for details, see:
+        # https://github.com/skyfielders/astronomy-notebooks
+        # and look for the notebook "error-in-timescale-ut1.ipynb".
+        tt = ut1 + delta_t_approx / DAY_S
+        t = Time(self, tt)
+        t.ut1 = ut1
+        return t
+
     def from_astropy(self, t):
         """Return a Skyfield time corresponding to the AstroPy time `t`."""
         return self.tt(jd=t.tt.jd)