Skip to content

Commit

Permalink
Add official routine for parsing Table S15
Browse files Browse the repository at this point in the history
  • Loading branch information
@brandon-rhodes
brandon-rhodes committed Mar 29, 2021
1 parent c41e7c6 commit ded5f52
Show file tree
Hide file tree
Showing 4 changed files with 145 additions and 131 deletions.
91 changes: 0 additions & 91 deletions Table-S15-v18.txt
View file

This file was deleted.

89 changes: 89 additions & 0 deletions Table-S15.2020.txt
View file
@@ -0,0 +1,89 @@
-+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Table S15: The Polynomial Coefficients for DT -720.0 to 2019.0 v. 2020
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Row Years Polynomial Coefficients
i K_i K_{i+1} a_0 a_1 a_2 a_3
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1 -720.0 -100.0 20371.848 -9999.586 776.247 409.160
2 -100.0 400.0 11557.668 -5822.270 1303.151 -503.433
3 400.0 1000.0 6535.116 -5671.519 -298.291 1085.087
4 1000.0 1150.0 1650.393 -753.210 184.811 -25.346
5 1150.0 1300.0 1056.647 -459.628 108.771 -24.641
6 1300.0 1500.0 681.149 -421.345 61.953 -29.414
7 1500.0 1600.0 292.343 -192.841 -6.572 16.197
8 1600.0 1650.0 109.127 -78.697 10.505 3.018
9 1650.0 1720.0 43.952 -68.089 38.333 -2.127
10 1720.0 1800.0 12.068 2.507 41.731 -37.939
11 1800.0 1810.0 18.367 -3.481 -1.126 1.918
12 1810.0 1820.0 15.678 0.021 4.629 -3.812
13 1820.0 1830.0 16.516 -2.157 -6.806 3.250
14 1830.0 1840.0 10.804 -6.018 2.944 -0.096
15 1840.0 1850.0 7.634 -0.416 2.658 -0.539
16 1850.0 1855.0 9.338 1.642 0.261 -0.883
17 1855.0 1860.0 10.357 -0.486 -2.389 1.558
18 1860.0 1865.0 9.040 -0.591 2.284 -2.477
19 1865.0 1870.0 8.255 -3.456 -5.148 2.720
20 1870.0 1875.0 2.371 -5.593 3.011 -0.914
21 1875.0 1880.0 -1.126 -2.314 0.269 -0.039
22 1880.0 1885.0 -3.210 -1.893 0.152 0.563
23 1885.0 1890.0 -4.388 0.101 1.842 -1.438
24 1890.0 1895.0 -3.884 -0.531 -2.474 1.871
25 1895.0 1900.0 -5.017 0.134 3.138 -0.232
26 1900.0 1905.0 -1.977 5.715 2.443 -1.257
27 1905.0 1910.0 4.923 6.828 -1.329 0.720
28 1910.0 1915.0 11.142 6.330 0.831 -0.825
29 1915.0 1920.0 17.479 5.518 -1.643 0.262
30 1920.0 1925.0 21.617 3.020 -0.856 0.008
31 1925.0 1930.0 23.789 1.333 -0.831 0.127
32 1930.0 1935.0 24.418 0.052 -0.449 0.142
33 1935.0 1940.0 24.164 -0.419 -0.022 0.702
34 1940.0 1945.0 24.426 1.645 2.086 -1.106
35 1945.0 1950.0 27.050 2.499 -1.232 0.614
36 1950.0 1953.0 28.932 1.127 0.220 -0.277
37 1953.0 1956.0 30.002 0.737 -0.610 0.631
38 1956.0 1959.0 30.760 1.409 1.282 -0.799
39 1959.0 1962.0 32.652 1.577 -1.115 0.507
40 1962.0 1965.0 33.621 0.868 0.406 0.199
41 1965.0 1968.0 35.093 2.275 1.002 -0.414
42 1968.0 1971.0 37.956 3.035 -0.242 0.202
43 1971.0 1974.0 40.951 3.157 0.364 -0.229
44 1974.0 1977.0 44.244 3.199 -0.323 0.172
45 1977.0 1980.0 47.291 3.069 0.193 -0.192
46 1980.0 1983.0 50.361 2.878 -0.384 0.081
47 1983.0 1986.0 52.936 2.354 -0.140 -0.165
48 1986.0 1989.0 54.984 1.577 -0.637 0.448
49 1989.0 1992.0 56.373 1.648 0.708 -0.276
50 1992.0 1995.0 58.453 2.235 -0.121 0.110
51 1995.0 1998.0 60.678 2.324 0.210 -0.313
52 1998.0 2001.0 62.898 1.804 -0.729 0.109
53 2001.0 2004.0 64.083 0.674 -0.402 0.199
54 2004.0 2007.0 64.553 0.466 0.194 -0.017
55 2007.0 2010.0 65.197 0.804 0.144 -0.084
56 2010.0 2013.0 66.061 0.839 -0.109 0.128
57 2013.0 2016.0 66.920 1.007 0.277 -0.095
58 2016.0 2019.0 68.109 1.277 -0.007 -0.139
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
The above table of polynomial coefficients enables evaluation of DT in
seconds (s) and its derivative (the length of day lod) in milliseconds
(ms) for any epoch between $-720$ and 2019. It is not valid outside the
specified range of years.

For the year and fraction Y, extract the coefficients a_0, a_1, a_2,
a_3 from row i, where K_i <= Y <= K_{i+1} and form
t = (Y - K_i)/(K_{i+1} - K_i), where 0 <= t < 1, and
thus calculate
DT = a_0 + a_1 t + a_2 t^2 + a_3 t^3 seconds
lod = (a_1 + 2 a_2 t + 3 a_3 t^2) / (K_{i+1}-K_i) / 0.36525 ms
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
These coefficients reproduce the spline approximation discussed by
L.V. Morrison, F.R. Stephenson, C.Y. Hohenkerk and M. Zawilski, in
their latest paper entitled ``Addendum 2020 to `Measurement of the
Earth's Rotation: 720 BC to AD 2015'' published in the Royal Society
Proceedings A, 478, 2021, see https://doi.org/10.1098/rspa.2020.0776.

Details of the original analysis is published in Royal Society
Proceedings A, 472, 2016, at https://doi.org/10.1098/rspa.2016.0404

All the data is also available from the website of HM Nautical Almanac
Office at http://astro.ukho.gov.uk/nao/lvm/.
-+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78 changes: 38 additions & 40 deletions design/delta_t.py
View file
Expand Up @@ -14,8 +14,8 @@
TODO:
[ ] Build reader for long-term file.
[ ] Build polynomial interpolation function.
[X] Build reader for long-term file.
[X] Build polynomial interpolation function.
[ ] Build translator between the two.
[ ] Build translator between finals2000A and interpolator.
[ ] Build combiner.
Expand All @@ -30,6 +30,7 @@
from skyfield import functions, timelib
from skyfield.api import load, Loader
from skyfield.data import iers
from skyfield.data.earth_orientation import parse_S15_table

inf = float('inf')

Expand All @@ -40,9 +41,11 @@ def main(argv):
#try_out_different_interpolation_techniques()

def compare_splines_to_finals2000_error_bars():
url = 'http://astro.ukho.gov.uk/nao/lvm/Table-S15-v18.txt'
#url = 'http://astro.ukho.gov.uk/nao/lvm/Table-S15-v18.txt'
url = 'http://astro.ukho.gov.uk/nao/lvm/Table-S15.2020.txt'
with load.open(url) as f:
columns = load_table_S15(f)
names, columns = parse_S15_table(f)

i, start_year, end_year, a0, a1, a2, a3 = columns

f = load.open('finals2000A.all')
Expand All @@ -51,32 +54,40 @@ def compare_splines_to_finals2000_error_bars():
iers._build_timescale_arrays(mjd_utc, dut1)
)

print(delta_t_recent.shape)
print(i.shape)
print('Size of IERS table:', delta_t_recent.shape)
print('Number of splines:', i.shape)

#year = [-720, 400, 700]
#year = np.arange(-720, 2010)
#year = np.arange(1800, 2010)
#year = np.arange(1980, 2010)
year = np.arange(1980, 2010, 0.1)
s15_curve = interpolate(year, start_year, a3, a2, a1, a0)
interpolate = prep_interpolate(start_year, a3, a2, a1, a0)
s15_curve = interpolate(year)

finals_tt, finals_delta_t = delta_t_recent
ts = load.timescale()
t = ts.utc(year)
tt = t.tt
print(tt)
print(finals_tt)
finals_curve = interpolate(
tt, finals_tt,

interpolate = prep_interpolate(
finals_tt,
finals_delta_t[1:] - finals_delta_t[:-1],
finals_delta_t[:-1],
)

print(finals_curve)
T0 = time()
finals_curve = interpolate(tt)
print(time() - T0, 's for interpolate()')

T0 = time()
finals_curve2 = np.interp(tt, finals_tt, finals_delta_t)
print(time() - T0, 's for interp()')

assert (finals_curve == finals_curve2).all()

diff = max(abs(s15_curve - finals_curve))
print('Max difference (s):', diff, diff * 15)
print('Max difference (seconds, arcseconds):', diff, diff * 15)

if 1:
fig, (ax, ax2) = plt.subplots(2, 1)
Expand All @@ -90,24 +101,20 @@ def compare_splines_to_finals2000_error_bars():
#plt.legend()
fig.savefig('tmp.png')

def interpolate(t, t_barriers, c, *coefficients):
def prep_interpolate(t_barriers, c0, *coefficients):
i_array = np.arange(len(t_barriers))
i = np.interp(t, t_barriers, i_array)
print(i)
# x = i % 1.0
# print(x)
i = i.astype(int)
print(i)
left = t_barriers[i]
right = t_barriers[i+1]
print(left, right)
little_t = (t - left) / (right - left)
print(little_t)
value = c[i]
for c in coefficients:
value *= little_t
value += c[i]
return value
def interpolate(t):
i = np.interp(t, t_barriers, i_array)
i = i.astype(int)
left = t_barriers[i]
right = t_barriers[i+1]
little_t = (t - left) / (right - left)
value = c0[i]
for c in coefficients:
value *= little_t
value += c[i]
return value
return interpolate

def try_adjusting_spline():
k0, k1 = -720.0, -100.0
Expand Down Expand Up @@ -210,7 +217,7 @@ def try_solving_spline():
def try_out_different_interpolation_techniques():
url = 'http://astro.ukho.gov.uk/nao/lvm/Table-S15-v18.txt'
with load.open(url) as f:
columns = load_table_S15(f)
names, columns = parse_S15_table(f)
i, start_year, end_year, a0, a1, a2, a3 = columns
report = []
print('Table start and end years:', start_year[0], end_year[-1])
Expand Down Expand Up @@ -475,14 +482,6 @@ def to_tt(year):

is_2050_to_2150_polynomial_worth_it()

def load_table_S15(f):
# http://astro.ukho.gov.uk/nao/lvm/Table-S15-v18.txt
content = f.read()
banner = b'- ' * 36 + b'-\n'
sections = content.split(banner)
table = np.loadtxt(sections[2].splitlines())
return table.T

def is_2050_to_2150_polynomial_worth_it():
# How different is the polynomial at:
# https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
Expand Down Expand Up @@ -565,4 +564,3 @@ def draw_plot_comparing_USNO_and_IERS_data():

if __name__ == '__main__':
main(sys.argv[1:])

18 changes: 18 additions & 0 deletions skyfield/data/earth_orientation.py
View file
@@ -1,3 +1,4 @@
# -*- coding: utf-8 -*-
"""Routines to download Earth orientation data."""

import numpy as np
Expand All @@ -13,6 +14,23 @@ def morrison_and_stephenson_2004_table():
df = tables[0]
return pd.DataFrame({'year': df[0], 'delta_t': df[1]})

def parse_S15_table(f):
"""Parse polynomial coefficients from Table S15.
The table is available at the website of Her Majesty's Nautical
Almanac Office, from the paper “Measurement of the Earth's Rotation:
720 BC to AD 2015” by L.V. Morrison, F.R. Stephenson, C.Y. Hohenkerk
and M. Zawilski 2021.
"""
# http://astro.ukho.gov.uk/nao/lvm/Table-S15.2020.txt
content = f.read()
banner = b'- ' * 36 + b'-\n'
sections = content.split(banner)
names = sections[1].splitlines()[-1].decode('utf-8').split()
table = np.loadtxt(sections[2].splitlines())
return names, table.T

def main():
thisdir = os.path.dirname(__file__)

Expand Down

0 comments on commit ded5f52

Please sign in to comment.