Repository: hydrogeoscience/pygtide (language: python)

Issue:
LOD tide is unsteady
The LOD tide shows unsteadiness. This results from the code interpolating data from etpolut1.dat and then differentiating it.
Is there a solution that makes the results steady?


Discussion / hints:
Can you provide the code snippet with which you generated the LOD and pole tidal contribution?
@hydrogeoscience I can work on this issue.
This is done in lines 2976-2999 in Fortran

Fix patch:
diff --git a/src/etpred.f90 b/src/etpred.f90
index dc1e275..df5bf38 100755
--- a/src/etpred.f90
+++ b/src/etpred.f90
@@ -658,7 +658,7 @@ SUBROUTINE PREDICT(ARGS)
       ITMIN=0
       ITSEC=0
 !-GCR fixed conversion
-!-GCR fixed total number of samples calculation 
+!-GCR fixed total number of samples calculation
       ! NDAT=NINT(DBLE(ISPANH)/DDTH)
       CALL ETJULN(IUN16,ITY,ITM,ITD,DTH,DJULD0)
 !#######################################################################
@@ -2808,6 +2808,67 @@ SUBROUTINE ETPHAS(IUN16,IPRINT,IMODEL,DLON,DJULD)
 END SUBROUTINE
 
 
+SUBROUTINE INTERP_QUADRATIC(DY1, DY2, DY3, DT, &
+        DOUT)
+!#######################################################################
+!     Original quadratic interpolation routine by Tom Eulenfeld
+!#######################################################################
+    IMPLICIT REAL(8) (D)
+
+    DOUT = DY2 + (DY3-DY1)*0.5D0*DT + (DY1-2.D0*DY2+DY3)*0.5D0 * DT * DT
+    ! difference is then given by
+    !DOUT = (DY3-DY1)*0.5D0 + (DY1-2.D0*DY2+DY3)* DT
+
+END SUBROUTINE
+
+
+SUBROUTINE INTERP(Y0, Y1, Y2, Y3, X0, DT, &
+        DOUT)
+!#######################################################################
+!     Cubic interpolation routine by Tom Eulenfeld
+!     Cubic spline interpolation of values between Y1 and Y2
+!     For cubic spline we need 4 constraints
+!     * curve goes through Y1
+!     * curve goes through Y2
+!     * derivative at Y1 given by differences Y2-Y0
+!     * derivative at Y2 given by differences Y3-Y1
+!#######################################################################
+    IMPLICIT REAL(8) (D)
+    REAL (8) x0, dx, y0, y1, y2, y3, a, b, c, d
+    dx = 1.0d0
+    a = 1.0d0*(2.0d0*dx**3*y0 - 3.0d0*dx**3*y1 + 3.0d0*dx**3*y2 - 1.0d0*dx** &
+        3*y3 + 4.0d0*dx**2*x0*y0 - 9.5d0*dx**2*x0*y1 + 8.0d0*dx**2*x0*y2 &
+        - 2.5d0*dx**2*x0*y3 + 2.5d0*dx*x0**2*y0 - 7.0d0*dx*x0**2*y1 + &
+        6.5d0*dx*x0**2*y2 - 2.0d0*dx*x0**2*y3 + 0.5d0*x0**3*y0 - 1.5d0*x0 &
+        **3*y1 + 1.5d0*x0**3*y2 - 0.5d0*x0**3*y3)/dx**3
+    b = 1.0d0*(-4.0d0*dx**2*y0 + 9.5d0*dx**2*y1 - 8.0d0*dx**2*y2 + 2.5d0*dx &
+        **2*y3 - 5.0d0*dx*x0*y0 + 14.0d0*dx*x0*y1 - 13.0d0*dx*x0*y2 + &
+        4.0d0*dx*x0*y3 - 1.5d0*x0**2*y0 + 4.5d0*x0**2*y1 - 4.5d0*x0**2*y2 &
+        + 1.5d0*x0**2*y3)/dx**3
+    c = 1.0d0*(2.5d0*dx*y0 - 7.0d0*dx*y1 + 6.5d0*dx*y2 - 2.0d0*dx*y3 + 1.5d0 &
+        *x0*y0 - 4.5d0*x0*y1 + 4.5d0*x0*y2 - 1.5d0*x0*y3)/dx**3
+    d = 1.0d0*(-0.5d0*y0 + 1.5d0*y1 - 1.5d0*y2 + 0.5d0*y3)/dx**3
+    DOUT = a + b * DT + c * DT * DT + d * DT * DT * DT
+END SUBROUTINE
+
+
+
+SUBROUTINE INTERP_DIF(DY0, DY1, DY2, DY3, DX0, DT, &
+        DOUT)
+!#######################################################################
+!     Differentiate and cubic interpolation routine by Tom Eulenfeld
+!#######################################################################
+    IMPLICIT REAL(8) (D)
+    REAL(8) Y0, Y1, Y2, Y3
+    Y1 = 0.5d0 * (DY2 - DY0)
+    Y2 = 0.5d0 * (DY3 - DY1)
+    Y0 = Y2 - 2.0d0*(DY2 - 2.0d0*DY1 + DY0)
+    Y3 = 2.0d0*(DY3 - 2.0d0*DY2 + DY1) + Y1
+    CALL INTERP(Y0, Y1, Y2, Y3, DX0, DT, DOUT)
+END SUBROUTINE
+
+
+
 SUBROUTINE ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, &
         DCLON,DSLON,DPOLX,DPOLY,DUT1,DTAI,DLOD,DGPOL,DGPOLP,DGLOD,NERR)
 !#######################################################################
@@ -2938,14 +2999,14 @@ SUBROUTINE ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, &
 !     DT is time difference referring to central sample point in days:
 !#######################################################################
       DT=DMODJD-DBLE(IMJD)
-      DT2=DT*DT
       IREC=IMJD-IFIRST+2
-      IF(IREC.LT.2) THEN
+      IF(IREC.LT.3) THEN
         DPOLX=0.D0
         DPOLY=0.D0
         DUT1 =0.D0
         DTAI =0.D0
         DLOD =0.D0
+        DGLOD =0.D0
         DGPOL=0.D0
         DGPOLP=0.D0
         NERR=1
@@ -2957,44 +3018,39 @@ SUBROUTINE ETPOLC(IUN16,IUN30,IUN31,IPRINT,DJULD,DCLAT,DSLAT, &
 !#######################################################################
         READ(IUN31,REC=ILAST)   DPOLX,DPOLY,DUT1,DTAI
         DLOD =0.D0
+        DGLOD =0.D0
         DGPOL=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLX*DCLON-DPOLY*DSLON)*DRAD/3600.D0*1.D9
         DGPOLP=0.D0
         NERR=1
         RETURN
       ENDIF
-      IF(IMJD.EQ.IMJDO) GOTO 1100
-	  !GCR fix bug here (which likely has no influence):
-	  ! READ(IUN31,REC=IREC-1) DPOLX1,DPOLY1,DUT12,DTAI1
       READ(IUN31,REC=IREC-1) DPOLX1,DPOLY1,DUT11,DTAI1
       READ(IUN31,REC=IREC)   DPOLX2,DPOLY2,DUT12,DTAI2
       READ(IUN31,REC=IREC+1) DPOLX3,DPOLY3,DUT13,DTAI3
-      IMJDO=IMJD
+      IF(IREC.EQ.ILAST-1) THEN
+        DPOLX4=DPOLX3
+        DPOLY4=DPOLY3
+        DTAI4=DTAI3
+      ELSE
+        READ(IUN31,REC=IREC+2) DPOLX4,DPOLY4,DUT14,DTAI4
+      ENDIF
 !#######################################################################
-!     Quadratic interpolation for pole coordinates and DTAI:
+!     Cubic interpolation for pole coordinates and DTAI:
 !     Linear interpolation for DUT1:
-!#######################################################################
- 1100 DPOLXA0=DPOLX2
-      DPOLXA1=(DPOLX3-DPOLX1)*0.5D0
-      DPOLXA2=(DPOLX1-2.D0*DPOLX2+DPOLX3)*0.5D0
-      DPOLYA0=DPOLY2
-      DPOLYA1=(DPOLY3-DPOLY1)*0.5D0
-      DPOLYA2=(DPOLY1-2.D0*DPOLY2+DPOLY3)*0.5D0
-      DTAIA0=DTAI2
-      DTAIA1=(DTAI3-DTAI1)*0.5D0
-      DTAIA2=(DTAI1-2.D0*DTAI2+DTAI3)*0.5D0
-      DUT10=DUT12
-      DDUT1=DUT13-DUT12
+!     TE: use cubic interpolation instead of quadratic
+!#######################################################################
+      CALL INTERP(DTAI1, DTAI2, DTAI3, DTAI4, -1.0d0, DT, DTAI)
+      CALL INTERP(DPOLX1, DPOLX2, DPOLX3, DPOLX4, -1.0d0, DT, DPOLX)
+      CALL INTERP(DPOLY1, DPOLY2, DPOLY3, DPOLY4, -1.0d0, DT, DPOLY)
+      CALL INTERP_DIF(DTAI1, DTAI2, DTAI3, DTAI4, -1.0d0, DT, DLOD)
+      CALL INTERP_DIF(DPOLX1, DPOLX2, DPOLX3, DPOLX4, -1.0d0, DT, DPOLXP)
+      CALL INTERP_DIF(DPOLY1, DPOLY2, DPOLY3, DPOLY4, -1.0d0, DT, DPOLYP)
+      DDUT1 = DUT13-DUT12
       IF(DDUT1.GT. 0.9D0) DDUT1=DDUT1-1.D0
       IF(DDUT1.LT.-0.9D0) DDUT1=DDUT1+1.D0
-      DLOD = DTAIA1+2.D0*DTAIA2*DT
+      DUT1 = DDUT1 * DT + DUT12
       DGLOD=2.D0*DLOD*DOM**2*DA*DCLAT*DCLAT*1.D9/86400.D0
-      DPOLX=DPOLXA0+DT*DPOLXA1+DT2*DPOLXA2
-      DPOLY=DPOLYA0+DT*DPOLYA1+DT2*DPOLYA2
-      DUT1 =DUT10  +DT*DDUT1
-      DTAI =DTAIA0 +DT*DTAIA1 +DT2*DTAIA2
       DGPOL=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLX*DCLON-DPOLY*DSLON)*DRAD/3600.D0*1.D9
-      DPOLXP=DPOLXA1+2.D0*DPOLXA2*DT
-      DPOLYP=DPOLYA1+2.D0*DPOLYA2*DT
       DGPOLP=DOM**2*DA*2.D0*DCLAT*DSLAT*(DPOLXP*DCLON-DPOLYP*DSLON)*DRAD/3600.D0*1.D9
       RETURN
 !#######################################################################


Test patch:
diff --git a/pygtide/tests.py b/pygtide/tests.py
index 2c21b08..6ade139 100755
--- a/pygtide/tests.py
+++ b/pygtide/tests.py
@@ -13,19 +13,19 @@ def test(msg=False):
     # test against a proven array
     args = (-20.82071, -70.15288, 830.0, '2017-01-01', 10, 600)
     series = predict_series(*args, statazimut=90, tidalcompo=0)
-    expected = np.array([ 607.605208,  630.191908,  646.733272,  657.225859,  661.710055,
-            660.269363,  653.029376,  640.156453,  621.856097,  598.371065,
-            569.979217,  536.991122,  499.747456,  458.616196,  413.989642,
-            366.28129 ,  315.922583,  263.35956 ,  209.049431,  153.457102,
-             97.051682,   40.30299 ,  -16.321912,  -72.362132, -127.366285,
-           -180.895731, -232.527698, -281.858255, -328.505135, -372.110372,
-           -412.342748, -448.900029, -481.510979, -509.937138, -533.974358,
-           -553.454086, -568.244386, -578.250699, -583.41634 , -583.72272 ,
-           -579.189308, -569.873325, -555.869181, -537.30765 , -514.354806,
-           -487.210711, -456.107887, -421.309551, -383.107667, -341.820786,
-           -297.791728, -251.38509 , -202.984615, -152.990441, -101.816229,
-            -49.886214,    2.367819,   54.509613,  106.101539,  156.707715,
-            205.89711 ])
+    expected = np.array([ 607.605208,  630.191906,  646.733265,  657.225843,  661.710028,
+        660.269321,  653.029317,  640.156373,  621.855993,  598.370935,
+        569.979057,  536.99093 ,  499.74723 ,  458.615933,  413.989338,
+        366.280944,  315.922193,  263.359123,  209.048945,  153.456565,
+        97.051092,   40.302345,  -16.322614,  -72.362893, -127.367107,
+    -180.896615, -232.528646, -281.859268, -328.506215, -372.11152 ,
+    -412.343966, -448.901318, -481.51234 , -509.938572, -533.975866,
+    -553.455669, -568.246045, -578.252435, -583.418153, -583.724611,
+    -579.191278, -569.875375, -555.87131 , -537.309859, -514.357095,
+    -487.213081, -456.110337, -421.312083, -383.110279, -341.823479,
+    -297.794502, -251.387943, -202.987549, -152.993454, -101.819321,
+        -49.889385,    2.36457 ,   54.506287,  106.098136,  156.704236,
+            205.893556])
     np.testing.assert_almost_equal(series, expected, 5)
 
     args = (-20.82071, -70.15288, 830.0, '2020-01-01', 29.5 * 24, 600)
@@ -45,6 +45,29 @@ def test(msg=False):
         plot_series(*args, statazimut=90, tidalcompo=8, show=False)
         if msg:
             plt.show()
-    
+
     print('------------------------------------')
     print('Successfully finished PyGTide tests!')
+
+
+def plot_pole_and_lod_tide():
+    try:
+        import matplotlib.pyplot as plt
+    except ImportError:
+        return False
+    ploty = ['Signal [nm/s**2]', 'Pole tide [nm/s**2]', 'LOD tide [nm/s**2]']
+    dates = ['1961-12-24', '2016-12-24', '2023-06-01']
+    for date in dates:
+        args = (20, 0, 0, date, 24*30, 600)
+        tab = predict_table(*args)
+        tab.plot(x='UTC', y=ploty, sharex=True, subplots=True)
+        plt.savefig(f'tides_{date}.png')
+    args = (-20.82071, -70.15288, 830.0, '2017-01-01', 10, 600)
+    tab = predict_table(*args, statazimut=90)
+    tab.plot(x='UTC', y=ploty, sharex=True, subplots=True)
+    plt.savefig('tides_2017-01-01_10h.png')
+    plt.show()
+
+
+
+#_plot_tides()