fix: use numpy interp for 2015 convention mean pole values

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'):