Switch to a more up-to-date data source for ∆T (and thus DUT1)

[guillaumechereau]

I am trying to compute dut1 using skyfield version 1.27:

import skyfield.api as sf
ts = sf.load.timescale()
t = ts.now()
print(t.dut1)

I get the value: -1.048504026495408 (I assume in seconds)

While I was expecting something closer to -0.176 seconds (source: http://hpiers.obspm.fr/eop-pc/index.php). The error seems quite big.

Is that an error in the delta_t computation, or maybe the leap seconds table? Or am I missing something else?

[brandon-rhodes]

Thanks for the question! The issue seems to be that Skyfield is using files that have fallen out of date, which it gets from here:

ftp://cddis.nasa.gov/products/iers/

Their deltat.data stopped being updated back in February, and ends with:

...
 2019 11  1  69.3432
 2019 12  1  69.3540
 2020  1  1  69.3612
 2020  2  1  69.3752

Once that file ends, Skyfield tries to use their "predictions" file for future dates, but we can see that as of the moment when they stopped updating those files, the predictions were already more than a year out of date — and predicted a much faster change than has actually occurred. Their deltat.preds starts with:

   58484.000  2019.00   69.34      -0.152       0.117       
   58575.000  2019.25   69.48      -0.295       0.162      
   58666.000  2019.50   69.62      -0.440       0.215    
   58758.000  2019.75   69.71      -0.527       0.273   
   58849.000  2020.00   69.87                   0.335       
   58940.000  2020.25   70.03                   0.399      
   59032.000  2020.50   70.16                   0.465       
   59123.000  2020.75   70.24                   0.532      
...

As you can see, they were predicting that ∆T would be at 70.2 by now, whereas in fact it's still down around 69.3 at this point. Which leads to the difference DUT1 you are looking at being off by exactly the same amount (nearly a full second!).

I suppose I should switch to Bulletin B for these values in the future:

ftp://hpiers.obspm.fr/iers/bul/bulb_new/

The inconvenience is that instead of being a single large file with the history going all the way back, it looks like I'll have to download and combine dozens of individual Bulletin B's in order to replace the single deltat.data that I used to use?

I probably will not be able to undertake that large a task this week, and then will be away from email and keyboards for all of next week. But I will update this issue when I get the chance to make some progress! Thanks for pointing out that it's apparently time to abandon the old data source.

[guillaumechereau]

No problem at all!

I looked a bit at the available files, and this one looks interesting:
ftp://cddis.nasa.gov/products/iers/finals.all
(description here: ftp://cddis.nasa.gov/products/iers/readme.finals)

It correctly predicts the DUT1 at the current date, have data starting from 1973, and is supposed to be updated weekly. It doesn't seem to directly contain ∆T though, but this could be computed from DUT1 and the leap second I suppose.

[brandon-rhodes]

Yes, that does look interesting! It will be annoying to use, because instead of relating the Earth's motion to a uniform timescale, it instead relates it only to the jagged UTC timescale — while making no indication at all about where leap seconds are, which I will have to take from another table.

Well, actually, that's not true, I could make an excellent guess about where leap seconds are because it's at the places where DUT1 suddenly jumps by nearly 1 second.

In any case, even though the actual physical property ∆T will now become a derived value, I suppose it happens behind the scenes anyway for Skyfield users. When I tackle this, probably in mid-October, I'll try using finals.all as Skyfield's new file-of-record.

Although, wow, 3.3 MB, versus about 14 kB for the current file! I don't imagine I'll want to carry the literal file around inside Skyfield like I do the ∆T files, it would make Skyfield 4× bigger than it currently is. Maybe I can boil down the file to a series of 4-byte records retaining only DUT1, which would take up 71k — much larger than the current file but still manageable.

I'll have to test how it affects performance; it will take longer for interpolation to find the correct daily value in the new file versus the monthly values from the old file. Maybe I should synthetically create a monthly extract? It would be the first time I'd broken with the strict layout of a fundamental data source and tried excerpting one instead, though, essentially creating my own new Skyfield-specific data source I'd be maintaining.

Anyway: lots of issues around the new file, but I'll see if it's a possibility the next time I have extended time to work on Skyfield!

[guillaumechereau]

I understand that a ∆T source is much nicer than a DUT1. In fact in the same NASA directory there is deltat.data and deltat.pred.
I also looked at the polynomial fitting done by Stephenson Morrison Hohenker (from -720 to 2013), and the algo from nasa website:
https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html, that is a description of the algo from Espenak and Meeus (2006).

Here is a graph comparing those values:

If you want I can send you the python script I used to generate this graph.

Not sure why the file deltat.pred doesn't match deltat.data even for the 2020 value.

[brandon-rhodes]

@guillaumechereau — All right! I think I have finished converting Skyfield over to the IERS for ∆T. Would you be interested in trying the feature out before the next release? You would need to:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

And then:

1. EITHER run a script that uses load.timescale() without the builtin= option,
2. OR run a script that uses load.timescale(builtin=False) after removing the three old files deltat.data, deltat.preds, and Leap_Second.dat from the file where you keep Skyfield download files.

If you do get a chance to test it out, let me know if Skyfield’s recent predictions for DUT1 look more reasonable. Thanks!

[guillaumechereau]

Thanks!

I just quickly added the skyfield computed value to my previous graph of different ∆T models. It seems to work as expected from 2000 to 2020. Something a bit strange seems to happen after 2022 though: the value goes up much faster than any other models.

[brandon-rhodes]

@guillaumechereau Thanks for trying it out! It looks like the old USNO deltat.preds had predictions all the way through year 2027.75, while the new finals2000A.all stops at 2021/12/12. Let me double-check the IERS and see if maybe there is a second separate predictions file I should be pulling it?

Once Skyfield reaches the end of the data file predictions, I think it just draws a straight line out to where the big curve from Morrison and Stephenson (2004) will be in 1 century from the end of the predictions, to try to provide a smooth transition between the two instead of an instantaneous jump somewhere.

I will draw a few graphs myself, and also go look for the papers you mention.

[brandon-rhodes]

@guillaumechereau — I have reviewed the old 2006 work of Espenak and Meeus at:

https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html

— and they do have a more elegant way of connecting the present with the long-term ∆T curve. Note the difference between how Skyfield currently just draws a straight line, and how they instead use a pair of curves:

Code that produces the above plot: https://github.com/skyfielders/python-skyfield/blob/a7298cbfea8205449f699a2b71ece16b1b622143/design/delta_t.py

I'm tempted to reduce their polynomials to, say, a yearly table of ∆T values up to 2150 that could be easily appended to the existing ∆T interpolation table. That would provide a better link than Skyfield's current straight line!

[guillaumechereau]

Thanks for the info.

For stellarium web (my original motivation to get a better ∆T) I am thinking about using a cubic spline interpolation during the time interval where there are measured values (like the Stephenson Morrison Hohenker 2016, but with coefficients recomputed to include the new observation data), and then do a parabolic curve fitting for the distant future, making sure it connects to the cubic spline.

An other thing: with skyfield, dut1 also seems to be large for distant past times:

t = ts.utc(-10000)
print(t.dut1)
>> -446999.0630375447

I also expected the value to say small, even though I understand that the definition of UTC might not really be defined at those times.

[brandon-rhodes]

Stephenson Morrison Hohenker 2016

Oh — I just opened that paper in a browser with cookies turned on, and had not realized that the full text of the article would appear! I’ll read it and see whether Skyfield should instead simply adopt their new model.

Or whether instead of invisibly switching between models, Skyfield should give users more control? Hmm.

An other thing: with skyfield, dut1 also seems to be large for distant past times:

The Morrison and Stephenson (2004) parabola does assume that the drift between clocks gets steeper and steeper, eventually generating quite a fast drift between TT and UT1.

[guillaumechereau]

The Morrison and Stephenson (2004) parabola does assume that the drift between clocks gets steeper and steeper, eventually generating quite a fast drift between TT and UT1.

Between TT and UT1 yes, but shouldn't UTC always be close to UT1 by design? That is, dut1 always within plus or minus one second?

[brandon-rhodes]

Oh — good question, I had misunderstood what you are asking. The answer is that the standards bodies have never published a retroactive schedule of leap seconds to be applied before 1972. Therefore, UTC before 1972 is simply another uniform time scale like TT and TAI, but fixed forever through history at TAI + 10s. It falls out of sync with historic UT1 at the same rate as the other uniform time scales.

[brandon-rhodes]

Though I hope to continue improving Skyfield's ∆T handling for dates before and after the ends of the IERS table, I'm going to go ahead and close this issue since with the recent release of Skyfield 1.31 the new IERS-powered ∆T is now in user’s hands, and no issues with it have been reported yet.

Thanks again for opening the issue and for the productive discussion! (And for the nice diagrams.)