diff --git a/timescale/eop.py b/timescale/eop.py index 498a047..16ee8b9 100644 --- a/timescale/eop.py +++ b/timescale/eop.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" eop.py -Written by Tyler Sutterley (04/2024) +Written by Tyler Sutterley (07/2025) Utilities for maintaining and calculating Earth Orientation Parameters (EOP) PYTHON DEPENDENCIES: @@ -14,6 +14,8 @@ utilities.py: download and management utilities for syncing files UPDATE HISTORY: + Updated 07/2025: use numpy interp for 2015 convention mean pole values + add Desai et al. (2015) secular pole model Updated 04/2024: fixed annotations with multiple types Updated 04/2023: using pathlib to define and expand paths add wrapper function for interpolating daily EOP values @@ -304,7 +306,7 @@ def update_finals_file( return # IERS mean or secular pole conventions -_conventions = ('2003', '2010', '2015', '2018') +_conventions = ('2003', '2010', '2015', 'Desai', '2018') # read table of mean pole values, calculate angular coordinates at epoch def iers_mean_pole( input_epoch: np.ndarray, @@ -313,7 +315,7 @@ def iers_mean_pole( ): """ Calculates the angular coordinates of the IERS Conventional Mean Pole (CMP) - or IERS Secular Pole (2018 convention) :cite:p:`Petit:2010tp` + or IERS Secular Pole (2018 convention) :cite:p:`Petit:2010tp,Desai:2015jr` Parameters ---------- @@ -326,6 +328,7 @@ def iers_mean_pole( - ``'2003'`` - ``'2010'`` - ``'2015'`` + - ``'Desai'`` - ``'2018'`` input_file: str or pathlib.Path Full path to mean-pole.tab file provided by IERS @@ -351,17 +354,11 @@ def iers_mean_pole( # read mean pole file input_file = pathlib.Path(kwargs['file']).expanduser().absolute() table = np.loadtxt(input_file) - # reduce to 1971 to end date - ii, = np.nonzero(table[:, 0] >= 1971) + # Reduce table following 2015 conventions: + # 1. trim dates prior to 1971 + # 2. only keep rows falling on exact years + ii, = np.nonzero((table[:, 0] >= 1971) & ((table[:, 0] % 1) == 0.0)) table = np.copy(table[ii,:]) - # reduce to yearly values - jj, = np.nonzero((table[:, 0] % 1) == 0.0) - table = np.copy(table[jj,:]) - end_time = table[-1, 0] + 0.2 - # final shape of the table - nrows, *_ = np.shape(table) - else: - end_time = np.inf # allocate for output arrays x = np.full_like(input_epoch, kwargs['fill_value']) y = np.full_like(input_epoch, kwargs['fill_value']) @@ -383,25 +380,19 @@ def iers_mean_pole( y[t] = 0.358891 - 0.0006287*dx flag[t] = True # Conventional mean pole model in IERS Conventions 2015 - # must be below maximum valid date within file (e.g. 2015.2 for 2015) - elif (convention == '2015') and (epoch >= 1975) and (epoch < end_time): - # find epoch within mean pole table - i = 1 - j = nrows+1 - while (j > (i+1)): - k = (i+j)//2 - if (epoch < table[k, 0]): - j = k - else: - i = k - # calculate differential from point in table - dx = epoch - table[i, 0] - if (i == (nrows-1)): - x[t] = table[i, 1] + dx*(table[nrows-1, 1]-table[nrows-2, 1]) - y[t] = table[i, 2] + dx*(table[nrows-1, 1]-table[nrows-2, 2]) - else: - x[t] = table[i, 1] + dx*(table[i+1, 1]-table[i, 1]) - y[t] = table[i, 2] + dx*(table[i+1, 2]-table[i, 2]) + # epoch must be within the dates in the mean pole file + elif (convention == '2015') and (epoch >= 1975): + # interpolate using times in table + x[t] = np.interp(epoch, table[:,0], table[:,1], + left=kwargs['fill_value'], right=kwargs['fill_value']) + y[t] = np.interp(epoch, table[:,0], table[:,2], + left=kwargs['fill_value'], right=kwargs['fill_value']) + flag[t] = (x[t] != kwargs['fill_value']) + # Secular pole model in Desai et al. (2015) + elif (convention == 'Desai'): + # calculate secular pole using equation 10a/b of Desai (2015) + x[t] = 0.05097 + 0.00062*(epoch - 2000.0) + y[t] = 0.33449 + 0.00348*(epoch - 2000.0) flag[t] = True # Secular pole model in IERS Conventions 2018 elif (convention == '2018'):