Values of ∆T = TT – UT1 overestimated for year > 2022

[miccoli]

The estimated values of ∆T seem to me to be greatly exaggerated for years > 2022.

Running this example

import skyfield, skyfield.api

print(skyfield.__name__, skyfield.__version__)
ts = api.load.timescale()
for yr in range(2021,2025):
    print("{}: ∆T = {}s".format(yr, ts.J(yr).delta_t))

yields

2021: ∆T = 69.35s
2022: ∆T = 69.78s
2023: ∆T = 71.80s
2024: ∆T = 73.82s

and of course it is not very likely to have an increment of ∆T of approx 2s / year in the near future. (This would amount to 2 leap seconds in 2022 and 2 leap seconds in 2023.)

I produced also some graphs to better analyse the problem.

LOD = length of day (which is the derivative of ∆T)
LOD - 86400s = difference between 1 SI day and 1 UT1 day

matplotlib program

tx = ts.J(np.linspace(2015, 2025, 200))

fig, ax = plt.subplots()
ax.plot(tx.J, tx.delta_t, label="∆t")
ax.legend()
ax.set_ylabel("(s)")
ax.set_xlabel("J(TT)")

Trend of ∆T in near future:

matplotlib program

tx = ts.J(np.arange(2000, 2200, 1/12))
dlod = np.diff(tx.delta_t) / np.diff(tx.ut1)

fig, ax = plt.subplots()
ax.plot(tx.J[1:], dlod*1000, label="LOD - 86400s")
ax.legend()
ax.axhline(0, color="red", alpha=0.5)
ax.set_ylabel("(ms)")
ax.set_xlabel("J(TT)")

Trend of LOD in near future

Long term trend of LOD:

I read #452, and it seems that current behaviour arises from the desire to "quickly" recover the extrapolated long term trend of ∆T (Stephenson, F.R., Morrison,L.V. and Hohenkerk, C.Y. (2016)).

But if we look at the resources in http://astro.ukho.gov.uk/nao/lvm/ it seems to me that we can expect that al least up to year 2500 the LOD will be well below the linear regression LOD obtained in years -2000 to 2000 (and this means that ∆T will also be below its long term trend.) On the contrary the Skyfield routines assume a dramatic step increase of the LOD in the near future, so that the long term ∆T will be almost recovered in the next century.

[brandon-rhodes]

I agree! Skyfield should, specifically:

• Adopt the splines listed under the heading “Observational Data and Spline Polynomial Coefficients for ΔT - Version 2018” at the http://astro.ukho.gov.uk/nao/lvm/ page.
• Offer the choice of historic/future ∆T model as something users can configure, both so that someone wanting an alternate curve can use it, and also so that applications can choose to explicitly pin their choice of the 2018 model so that future Skyfield updates don’t “improve” their results without their intending it.
• Another configuration choice should be the “old behavior” of the current table and parabola based on Morrison and Stephenson (2004) parabola, so that folks can choose to keep getting the Skyfield results they're used to seeing if they want them.
• Users should be able to configure how Skyfield sutures together the different sources of ΔT. Unless a user wants an abrupt jump, there will have to be a transition (a) from the long-term parabola to the first spline, (b) from the splines to the first data point of finals2000A.all, (c) from the end of finals2000A.all and back to the splines, and then (d) from the final spline on to the long-term parabola. (Or do the 2018 splines already provide smooth transitions for (a) and (d)?) And, we will have to choose a default way to smoothly join them for users who aren't interested in choosing a mechanism of their own. How long should the transitions be? How should they behave? We'll have to choose a reasonable default.

It would be nice, while accomplishing all that, to not slow down Skyfield too measurably, even though spline interpolation will be more expensive than linear interpolation. If you would like to see my work on the problem so far, check out:

https://github.com/skyfielders/python-skyfield/blob/master/design/delta_t.py

You can see that I made some progress back in October, but then switched projects to polar motion — once I'd written the code to parse finals2000A.all, which also provides the polar motion offsets, I figured that I should go ahead and tackle a longstanding issue that involved polar motion.

Thanks for opening this issue, it provides a useful way to track further progress. And feel free to weigh in with any design thoughts of your own towards the new spline-interpolation feature!

[miccoli]

I will for sure look at your work on delta_t.py. Let me share just a few random thoughts.

If I got it right, the sources of ∆T in http://astro.ukho.gov.uk/nao/lvm/ are

1. The "average parabola" (linear LOD)
2. The extrapolated long term trend (linear LOD + sine term)
3. The "spline approximation" valid for -720 < year < 2016 (piecewise quadratic LOD)

(All these approximations are obtained by a least squares fit of the historical data available from -720 and forward.)

All these expression give raise to different ∆T and have a different associated error, so it makes little sense to me to stitch them together. The user should therefore choose the most appropriate solution, depending on the application and the range of years of interest.

If an "automagic" value is desired, the library should

• for point values, choose the best approx available
• for ranges, choose the approx that covers the entire range requested.

[miccoli]

I have some preliminary comparisons of 1., 2., 3., vs Skyfield ts.delta_t in terms of LOD. (LOD obtained analytically for 1., 2., 3., and by 2 point finite differences for Skyfield).

Excess LOD in ms, as a function of year:

Comparison in terms of ∆T will follow.

[miccoli]

I've computed also estimates of ∆T.

As what regards 2. (long term extrapolation), it is obtained by integration of the corresponding LOD expression. I've only skimmed through the paper, but by comparison with the figures in http://astro.ukho.gov.uk/nao/lvm/ it seems to me that the integration constant has to be determined by continuity with the spline approximation. Moreover it seems that this expression should only be integrated forward, and not backward in time. Nevertheless here are some comparisons.

The overall picture looks good:

However a zoom on the 1500–2200 period (omitting parabola 1) shows that Skiefield extrapolation is not satisfactory:

If we look at the 2000–2500 extrapolation this is the picture:

Finally the difference between all approximations and parabola 1. (hidden because very distracting)

Difference of ∆T formulas with respect to parabola 1. as baseline

Conclusion

What to do for values of ∆T in the past, I don't have a strong opinion now. But for extrapolation in the future, I have no doubt that the most sensible approach would be to integrate the expression 2. of the LOD starting from the last available point in the IERS file.

If you agree, I can provide a PR in this sense:

• no changes in the current approach for values in the past
• for values after the last point in finals.all, integrate expression 2. of the LOD.

[brandon-rhodes]

(I'm marking this as a "feature request" to remind myself that it will require new code, to supplement the current long-term ∆T with more polynomial segments.)

[brandon-rhodes]

@miccoli — I might have a bit of time in my schedule now to dig into this a bit further and more fully respond. At the http://astro.ukho.gov.uk/nao/lvm/ page, I now see:

• Addendum 2020 to `Measurement of the Earth's Rotation: 720 BC to AD 2015', published in 2021 February in the Royal Society's Proceedings A 478, by Morrison, L. V., Stephenson, F.R., Hohenkerk, C.Y. and M. Zawilski, M.

So it looks like there are updated polynomials to use for the past ∆T. It looks like for length of day from 2020 to 2500, they now integrate this expression:

lod = +1.72 t − 3.5 sin(2π(t+0.75)/14)

where t is the number of centuries from +1825.0. I kind of wish they had a spline instead, for symmetry with the way they handle the years earlier than that.

What approach were you thinking of — were you going to build a linear interpolation table out into the future with the integration? Or use the formula itself?

It's an interesting idea to put the integrated curve out at the end of the final observation in the IERS file. Do we really want someone's code that's, say, predicting solar eclipses to change the details of its predictions with every little wiggle of the most recent day's IERS number — which then changes over the weeks that follow as they do smoothing? Or should the integration constant "anchor" be one that wiggles around a bit less?

[brandon-rhodes]

Are my eyes deceiving me, or do Morrison, Stephenson, Hohenkerk, and Zawilski themselves choose to position their long-term curve independently of their final interpolated value, instead of integrating forward from it? It looks to me like their long-term curve does not meet the tip of their final spline:

(From http://astro.ukho.gov.uk/nao/lvm/ in the 2020–2500 section)

[miccoli]

This is the image captured from https://web.archive.org/web/20200218152225/http://astro.ukho.gov.uk/nao/lvm/ which I analysed for "guessing" the integration constant for the 2018 long term ∆T prediction.

The 2018 expression of the LOD was

lod = + 1.78 t - 4.0 sin(2πt/15)

Probably my assumption of continuity was wrong: looking very closely there is a hint that the last spline point (2016) is slightly above the extrapolated curve (2010–2200). But somehow in the 2010–2016 interval both the spline and the interpolated curve appear to be very close, so it is hard to guess how the integration constant was really determined. However the graph that I obtained by continuity with the last spline point was remarkably in good agreement with the 2018 Morrison et.al. graph.

But with the 2021 graph, there is no doubt that my assumption was wrong!

As what regards the new situation, I had no time to properly analyse the new data but here are some quick remarks.

• Discontinuities of ∆T should be avoided. But I agree that having a prediction of ∆T that is very sensible to the last IERS data points is not a good idea. Somehow high frequency oscillations has to dampened, or filtered out. From this point of view having a fixed long term prediction seems sensible, provided that a bridge function is devised, so that short term to medium term predictions are computed, which are continuous with both the last IERS data and the long term predictions.
• My original idea was to use the analytical expression of the long term ∆T, which is a parabola plus a cosine term.
• Probably a mixed approach could be used: long term expression of ∆T is analytical, short term to medium term prediction is a set of splines, whose coefficients are computed when new IERS data are loaded. (User configurable, so that old predictions can be kept unchanged, at the cost of a discontinuity in the ∆T.)

But for now, I have not yet elaborated on how to bridge between IERS data and long term prediction.

[brandon-rhodes]

Well, that took a few weeks longer than I had expected, but after investigating the options I have just released Skyfield 1.38 with what I hope you will find is a gentler ∆T curve over the next few decades.

Note that I opted not to experiment with the new sin() term, since it wasn’t present in the 2016 open-access paper, and since I haven't seen anything since their 2016 paper that gives the integrated long-term parabola to reveal their choice of integration constant — and also, the code comes out both simpler and faster. Let’s see if their next couple of papers also keep including the sin() term; if it winds up sticking around, and its parameters are fairly stable, then we could at least update the Skyfield docs to show people how to add it in, even if we didn't make it the default because of the extra computational expense.

Let me know if you see any problems, or have difficultly getting out a curve that looks like the one above!

[miccoli]

Kudos for the new code, which is open to a lot of possible user customisations!

First a general remark on the sin() term. After reading again the 2016 paper I agree that it makes no sense to introduce this term. From the 2016 paper, page 22:

The reality of the quasi-periodic behaviour of the lod should be treated guardedly.

and

The extrapolation of the lod beyond the limits of the dataset is dependent on the reality of the 1500 year oscillation. It is also assumed that tidal friction is essentially unchanged to within the boundaries of its uncertainty in the period −2000 to 2500. Both are somewhat conjectural and require further investigation in the future.

I initially thought that the sin term was desirable because of the remarkably smooth transition from the last spline to the extrapolated (parabola + sin) function. But it turns out this was only by chance! Also the idea to fix the integration constant by continuity with the last known point turned out to be quite simplistic. Assimilation of the last IERS data points should be done in a more comprehensive way (by a Kalman filter?).

So I agree that the skyfield v1.38 implementation is a very sensible one.

Just a few remarks on some minor details.

1. skyfield v1.38 is using the 2016 parabola (∆T = 32.5 * t**2 - 320, lod = 1.78 * t) with the 2020 splines of http://astro.ukho.gov.uk/nao/lvm/. However the linear term of the 2020 lod is 1.72 * t which when integrated becomes 31.4 * t**2 + X where again X is an unknown integration constant. I assume this choice was made because there is no open access documentation for the 2020 parabola, but does it make sense to mix the 2016 long term parabola with the 2020 splines?
2. There is a small discontinuity of ∆T between the splines and the first point of the IERS tables (1973-01-02T00:00UTC).

example code

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import skyfield
import skyfield.api

ts = skyfield.api.load.timescale(builtin=True)
tx = ts.utc(
    year=1973, month=1, day=2, hour=np.linspace(-24, 24, 1000), minute=0, second=0
)

fig, ax = plt.subplots()
ax.plot(tx.utc_datetime(), tx.delta_t, label="Time.delta_t")
ax.plot(
    tx.utc_datetime(),
    ts.delta_t_function.long_term_function(tx.J),
    linestyle="dashed",
    label="long_term_function",
)
ax.grid()
ax.legend()
ax.set_xlabel("UTC")
ax.set_ylabel("∆T (s)")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%dT%H:%M"))
fig.autofmt_xdate()

[brandon-rhodes]

First a general remark on the sin() term…

The sin() term would have been fun to include! Though it’s tentative, it’s a very good match for the shape of their lod curve. I mainly avoided it because the long-term parabola can be very simply rolled into the Splines object, whereas applying the sin() would mean turning the 3 steps interpolate / detect-nans / compute-splines into 5 steps: interpolate / detect-nans / compute-splines / detect-nans-again / compute-parabola-plus-sine.

I assume this choice was made because there is no open access documentation for the 2020 parabola, but does it make sense to mix the 2016 long term parabola with the 2020 splines?

I can't think of any problems offhand, as they are both good sets of numbers based on very recent analysis, but if you come across any problems I'll definitely want to hear about them.

There is a small discontinuity of ∆T between the splines and the first point of the IERS tables (1973-01-02T00:00UTC).

Drat! I forgot about that one, and the graph I was using to look for discontinuities obscured that one since it was right up against the noise from the IERS daily values. As you can see above, I've committed a fix; if you'd like to try it out yourself ahead of the next Skyfield release, you can:

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