diff --git a/Table-S15-v18.txt b/Table-S15-v18.txt deleted file mode 100644 index 9e3521ff7..000000000 --- a/Table-S15-v18.txt +++ /dev/null @@ -1,91 +0,0 @@ --+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -Table S15: The Polynomial Coefficients for DT -720.0 to 2016.0 v. 2018 -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -Row Years Polynomial Coefficients - i K_i K_{i+1} a_0 a_1 a_2 a_3 -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1 -720.0 400.0 20483.158 -20530.285 10289.578 -3726.121 - 2 400.0 1000.0 6516.330 -5962.228 -255.072 1244.067 - 3 1000.0 1300.0 1543.097 -1370.086 869.282 -378.827 - 4 1300.0 1500.0 663.466 -512.002 -118.755 169.232 - 5 1500.0 1600.0 201.941 -120.908 97.235 -65.127 - 6 1600.0 1650.0 113.140 -60.910 -24.537 17.734 - 7 1650.0 1720.0 45.427 -79.494 56.185 -9.757 - 8 1720.0 1800.0 12.360 4.118 35.151 -33.363 - 9 1800.0 1810.0 18.267 -3.209 -1.015 1.752 - 10 1810.0 1820.0 15.795 0.017 4.240 -3.582 - 11 1820.0 1830.0 16.471 -2.247 -6.504 3.097 - 12 1830.0 1840.0 10.816 -5.965 2.787 -0.014 - 13 1840.0 1850.0 7.625 -0.432 2.745 -0.594 - 14 1850.0 1855.0 9.343 1.637 0.241 -0.870 - 15 1855.0 1860.0 10.352 -0.490 -2.368 1.547 - 16 1860.0 1865.0 9.040 -0.586 2.272 -2.472 - 17 1865.0 1870.0 8.254 -3.457 -5.143 2.717 - 18 1870.0 1875.0 2.371 -5.593 3.008 -0.913 - 19 1875.0 1880.0 -1.126 -2.315 0.270 -0.039 - 20 1880.0 1885.0 -3.210 -1.893 0.152 0.563 - 21 1885.0 1890.0 -4.388 0.101 1.842 -1.439 - 22 1890.0 1895.0 -3.884 -0.531 -2.474 1.871 - 23 1895.0 1900.0 -5.017 0.134 3.138 -0.232 - 24 1900.0 1905.0 -1.977 5.715 2.443 -1.257 - 25 1905.0 1910.0 4.923 6.828 -1.329 0.720 - 26 1910.0 1915.0 11.142 6.330 0.831 -0.825 - 27 1915.0 1920.0 17.479 5.518 -1.643 0.262 - 28 1920.0 1925.0 21.617 3.020 -0.856 0.008 - 29 1925.0 1930.0 23.789 1.333 -0.831 0.127 - 30 1930.0 1935.0 24.418 0.052 -0.449 0.142 - 31 1935.0 1940.0 24.164 -0.419 -0.022 0.702 - 32 1940.0 1945.0 24.426 1.645 2.086 -1.106 - 33 1945.0 1950.0 27.050 2.499 -1.232 0.614 - 34 1950.0 1953.0 28.932 1.127 0.220 -0.277 - 35 1953.0 1956.0 30.002 0.737 -0.610 0.631 - 36 1956.0 1959.0 30.760 1.409 1.282 -0.799 - 37 1959.0 1962.0 32.652 1.577 -1.115 0.507 - 38 1962.0 1965.0 33.621 0.868 0.406 0.199 - 39 1965.0 1968.0 35.093 2.275 1.002 -0.414 - 40 1968.0 1971.0 37.956 3.035 -0.242 0.202 - 41 1971.0 1974.0 40.951 3.157 0.364 -0.229 - 42 1974.0 1977.0 44.244 3.198 -0.323 0.172 - 43 1977.0 1980.0 47.291 3.069 0.193 -0.192 - 44 1980.0 1983.0 50.361 2.878 -0.384 0.081 - 45 1983.0 1986.0 52.936 2.354 -0.140 -0.166 - 46 1986.0 1989.0 54.984 1.577 -0.637 0.448 - 47 1989.0 1992.0 56.373 1.649 0.709 -0.277 - 48 1992.0 1995.0 58.453 2.235 -0.122 0.111 - 49 1995.0 1998.0 60.677 2.324 0.212 -0.315 - 50 1998.0 2001.0 62.899 1.804 -0.732 0.112 - 51 2001.0 2004.0 64.082 0.675 -0.396 0.193 - 52 2004.0 2007.0 64.555 0.463 0.184 -0.008 - 53 2007.0 2010.0 65.194 0.809 0.161 -0.101 - 54 2010.0 2013.0 66.063 0.828 -0.142 0.168 - 55 2013.0 2016.0 66.917 1.046 0.360 -0.282 -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -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 2016. 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 and C.Y. Hohenkerk, in ``Historical -Changes in UT/lod from Eclipses", Proceedings of the Journees des -Systemes de Reference et de la Rotation Terrestre, 2017, -https://web.ua.es/journees2017/ (in preparation). -This is an update of our Royal Society paper (2016)(see below) for the -inclusion of three hitherto excluded Chinese eclipses in 198 BC, AD 306 -and AD 616, and two critical Spanish observations in AD 1239 and -strengthened the case for the Chinese observation AD 1361. - -The details of the observations and their reduction are available in our -Royal Society paper and in the supplementary material: ``Measurement of -the Earth's Rotation: 720 BC to AD 2015'', Stephenson, F.R., Morrison, -L.V., and Hohenkerk, C.Y., published by Royal Society Proceedings A, -and available from their website at -http://rspa.royalsocietypublishing.org/lookup/doi/10.1098/rspa.2016.0404 --+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/Table-S15.2020.txt b/Table-S15.2020.txt new file mode 100644 index 000000000..3fa51132b --- /dev/null +++ b/Table-S15.2020.txt @@ -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/. +-+- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/design/delta_t.py b/design/delta_t.py index ccd646cf8..cf73b7e03 100755 --- a/design/delta_t.py +++ b/design/delta_t.py @@ -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. @@ -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') @@ -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') @@ -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) @@ -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 @@ -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]) @@ -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 @@ -565,4 +564,3 @@ def draw_plot_comparing_USNO_and_IERS_data(): if __name__ == '__main__': main(sys.argv[1:]) - diff --git a/skyfield/data/earth_orientation.py b/skyfield/data/earth_orientation.py index 9dd3ea8b3..d6cdeca1d 100644 --- a/skyfield/data/earth_orientation.py +++ b/skyfield/data/earth_orientation.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """Routines to download Earth orientation data.""" import numpy as np @@ -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__)