Prominently document that UT1, not UTC, is most appropriate for ancient dates

[jacob-pro]

I am getting very strange results for sunset, the older the date, it seems to be out by a few hours, (Expected value from Starry Night software)

Got 2020-06-01T16:05:01Z, expected 2020-06-01T16:05:00Z, diff mins 0, alt -0.83
Got 1800-06-01T16:03:42Z, expected 1800-06-01T16:04:00Z, diff mins 0, alt -0.83
Got 1000-06-01T16:31:04Z, expected 1000-06-01T16:05:00Z, diff mins 26, alt -0.83
Got 0200-06-01T18:23:57Z, expected 0200-06-01T16:00:00Z, diff mins 143, alt -0.83
Got -400-06-01T20:14:49Z, expected -400-06-01T15:56:00Z, diff mins 258, alt -0.83
Got -800-06-01T21:59:03Z, expected -800-06-01T15:53:00Z, diff mins 366, alt -0.83

To reproduce:

load = Loader(path.as_posix())
ts = load.timescale()
ts.julian_calendar_cutoff = GREGORIAN_START
eph = load('de431t.bsp')
babylon = Topos("32.55 N", "44.42 E")

def altitude_of_sun(t0):
    b = eph['earth'] + babylon
    apparent = b.at(t0).observe(eph['sun']).apparent()
    alt, az, distance = apparent.altaz()
    return alt

def sunset_for_day(t0, expected):
    t1 = ts.tt_jd(t0.tt + 1)
    t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(eph, babylon))
    for ti, yi in zip(t, y):
        if yi == 0:
            diff_mins = int(abs(ti.tt - expected.tt) * 60 * 24)
            print("Got {}, expected {}, diff mins {}, alt {:.2f}"
                  .format(ti.utc_iso(), expected.utc_iso(), diff_mins, altitude_of_sun(ti).degrees))

if __name__ == "__main__":
    sunset_for_day(ts.utc(2020, 6, 1), ts.utc(2020, 6, 1, 16, 5))
    sunset_for_day(ts.utc(1800, 6, 1), ts.utc(1800, 6, 1, 16, 4))
    sunset_for_day(ts.utc(1000, 6, 1), ts.utc(1000, 6, 1, 16, 5))
    sunset_for_day(ts.utc(200, 6, 1), ts.utc(200, 6, 1, 16, 0))
    sunset_for_day(ts.utc(-400, 6, 1), ts.utc(-400, 6, 1, 15, 56))
    sunset_for_day(ts.utc(-800, 6, 1), ts.utc(-800, 6, 1, 15, 53))

Am I doing something wrong? I'm guessing its something to do with the timescale? I'm on version 1.34

[brandon-rhodes]

(Merry Christmas!) The UTC timescale uses leap seconds to keep up with the Earth's actual rotation, but leap seconds are only specified from 1972 forwards — there are no leap seconds going backwards that would keep it in sync in past centuries with the length of the day, or at least with our retrospective guesses about how long the day was in past centuries. We currently use the sparse ancient records of eclipses to try to guess where the Earth was pointed during the day in ancient times, and those guesses guide the UT1 timescale that follows the Sun rather than an artificial scheme of leap seconds. I'd first try UT1 and see whether that gives you closer values to the ones you expect.

[jacob-pro]

Merry christmas! Thanks so much! I thought i had tried that but clearly not. Replacing all utc() with ut1() seems to work as expected

Got 1000-06-01T16:31:04Z, expected 1000-06-01T16:30:26Z, diff mins 0, alt -0.83
Got 0200-06-01T18:23:57Z, expected 0200-06-01T18:23:14Z, diff mins 0, alt -0.83
Got -400-06-01T20:14:49Z, expected -400-06-01T20:14:02Z, diff mins 0, alt -0.83
Got -800-06-01T21:59:03Z, expected -800-06-01T21:57:56Z, diff mins 1, alt -0.83

[brandon-rhodes]

Great, I'm glad it helped! If you don't mind, I'm going to re-open this with a different title to remind myself to highlight this difference more in the documentation.

[miccoli]

Let me add a small note. I'm new to this (truly wonderful, thanks!) package, and I was a little surprised by the utc_ methods.

For example look at

>>> from skyfield import api
>>> ts = api.load.timescale(builtin=True)
>>> ts.ut1(0).ut1_strftime()
'0000-01-01 00:00:00 UT1'
>>> ts.ut1(0).utc_jpl()
'B.C. 0001-Jan-01 02:55:37.8127 UT'

The utc_jpl() method produces a timestamp under these two assumptions:

1. UTC before 1972 ticks at 1 SI second frequency
2. 'UT' in the JPL timestamp refers to this "proleptic" UTC.

But 1. is not historically accurate: the term 'Coordinated Universal Time (UTC)' was already used in the sixties, referring to a radio broadcast time/frequency signal, ticking with a uniform rate (quartz oscillators, whose long term drift was controlled by caesium clocks) and that was kept within 50 ms to UT2 by 50ms step adjustments. This radio signal was first standardised in 1963 by the CCIR, and evolved in the current UTC, standardised by ITU-R in 1972, and linked to TAI.

As what regards 2., "UT" in the JPL time stamp, refers to the concept of universal time (UT), which is "mean solar time referred to the Greenwich meridian counted from midnight". (This because until 1925, Greenwich mean time (GMT) was used in astronomy for the mean solar time counted from noon, and a different name was needed for time counted from midnight.)

More precisely, UT in the JPL timestamp is defined in this document as

UT is Universal Time This can mean one of two non-uniform time-scales based on the rotation of the Earth. For this program, prior to 1962, UT means UT1. After 1962, UT means UTC or "Coordinated Universal Time". Future UTC leap-seconds are not known yet, so the closest known leap-second correction is used over future time-spans.

The time stamp produced by Skyfield (B.C. 0001-Jan-01 02:55:37.8127 UT) is therefore highly misleading for me, since it refers to a timescale that I'm unaware of.

[brandon-rhodes]

Thanks for your informative comment, @miccoli! I think it can result in several improvements.

1. I would not be opposed to improving Skyfield's definition of UTC back before 1972, if we could get hold of a table of those 50 ms adjustments. Do you know if a schedule of them exists? I was not aware of any use of the acronym "UTC" before 1970, and am very interested in any contemporary 1960s documentation for UTC that we could find.
2. I agree that the output of utc_jpl() is misleading — because I copy and pasted from the HORIZONS text output format without adapting to the actual calculation Skyfield is performing, it wrongly says UT instead of UTC. I will adjust the output to say explicitly UTC.
3. I believe that we should keep the behavior that methods with utc in the name return UTC, not any other UT.
4. But we should probably introduce a method ut_jpl(), and maybe matching methods ut_strftime() and ut_calendar() and so forth, that use UTC for recent decades but fall back gracefully to estimated UT1 for ancient dates and times. Currently users have to choose manually between when to ask for UTC and when to ask for UT1. It would be convenient to have a reasonable default built in to Skyfield that ut_…() methods use; it could be configurable on the timescale object exactly how the cutover is performed.
5. See #519 for discussion of improvements currently planned for Skyfield’s UT1 curves. That issue discusses how various curves are joined to form UT1, and how that could be configurable by users. Maybe the same mechanism could be available — with a reasonable built-in default — for users to configure how UTC and UT1 are joined together to form UT.

Do you think that would satisfy your use case?

[miccoli]

(Disclaimer: I'm not a professional astronomer, so my opinions here are not "authoritative")

1. Useful (academic) sources on the history of UTC are
   • The leap second: its history and possible future, R A Nelson et al 2001 Metrologia 38 p. 509
   • Mean Solar Time on the Meridian of Greenwich, DH Sadler, 1978 Q. Jl R. astr. Soc. 19 p. 290

According to these sources the first "version" of UTC was formalised by

• International Radio Consultative Committee (CCIR), Recommendation 374, Standard-Frequency and Time- Signal Emissions, Documents of the Xth Plenary Assem- bly, Geneva, Switzerland, 1963, Geneva, International Telecommunication Union, 1963, Vol. III, 193.

Quoting Nelson 2001:

The nautical almanacs of the UK and the USA were combined in 1957, beginning with the editions for 1960. In August 1959 it was also agreed to coordinate their time and frequency transmissions. Coordination began 1 January 1960. The participating observatories and laboratories were the USNO, RGO, NBS, NRL and NPL. Gradually other countries joined the system, which was entrusted to the BIH in 1961. In January 1965, the BIH decided to attach UTC to its atomic time A3 (which became TAI) by a mathematical relationship [84]. This was the origin of the link between TAI and UTC. The name “Coordinated Universal Time (UTC)” was approved by a resolution of IAU Commissions 4 and 31 at the 13th General Assembly in 1967 [85] .... Details of the UTC system were formalized by CCIR Study Group 7 in Geneva in 1962 and were adopted by the CCIR in its Recommendation 374 [86] of 1963.

Unfortunately I have no primary sources on the 1963-1972 UTC – UT1 difference. (I assume 50ms corrections up to 1962, then 100ms corrections to remain in sync with UT2, which is within 35ms of UT1, so the difference was smaller than the present DUT1). Authoritative data should be available at the former BIH (Bureau International de l'Heure). If I find some source I will link it here.

2. Good solution. Maybe the docs should state that for dates before 1972, (or 1960 if we find the relevant offsets) the results returned by Skyfield refer to a timescale that scale was not in use, and may be not accurate.

3. and 4. I agree. ut_jpl and ut_strftime would be very useful, since this time (with a UT1 to UTC transition) reflects the notion of clock-time that most people have for the civil use. I mean, one would set its clock to UT1 pre 1962, and to UTC (radio signals) afterwards. And projecting UT in the ancient past is still useful, although the ancient time scales were based on the motion of the apparent sun and not the mean sun. (May be a function for the equation of time could be a good addition to Skyfield... but this is a different subject.)

4. Yes: a sound default, but with the possibility of user configuration is the best option. Personally I would like also the option of the utc_ methods raising an exception if used when UTC was not yet defined.

[brandon-rhodes]

Personally I would like also the option of the utc_ methods raising an exception if used when UTC was not yet defined.

That’s an interesting idea! And it would be easy for the docs to show how. I am currently working to improve the internal leap-second tables, at which point the following proof-of-concept might stop working, but with the current code here’s how it could be implemented by an end user today:

import numpy as np
from skyfield.api import load

ts = load.timescale()
print(ts.leap_dates)
print(ts.leap_offsets)

ts.leap_offsets[0:2] = np.nan

t = ts.utc(1972, 1, 1)
print(t.utc_jpl())

t = ts.utc(1971, 12, 30)
print(t.utc_jpl())

That returns a ValueError for the second date when I try it.