refactor: make offsets to convert to MJD variables (#306)

tsutterley · Jul 9, 2024 · 3c8789f · 3c8789f
1 parent 747e7f4
commit 3c8789f
Show file tree

Hide file tree

Showing 4 changed files with 43 additions and 21 deletions.
diff --git a/pyTMD/astro.py b/pyTMD/astro.py
@@ -17,6 +17,7 @@
 
 UPDATE HISTORY:
     Updated 07/2024: made a wrapper function for normalizing angles
+        make number of days to convert days since an epoch to MJD variables
     Updated 04/2024: use wrapper to importlib for optional dependencies
     Updated 01/2024: refactored lunisolar ephemerides functions
     Updated 12/2023: refactored phase_angles function to doodson_arguments
@@ -63,6 +64,13 @@
 # default JPL Spacecraft and Planet ephemerides kernel
 _default_kernel = get_data_path(['data','de440s.bsp'])
 
+# number of days between the Julian day epoch and MJD
+_jd_mjd = 2400000.5
+# number of days between MJD and the J2000 epoch
+_mjd_j2000 = 51544.5
+# Julian century
+_century = 36525.0
+
 # PURPOSE: calculate the sum of a polynomial function of time
 def polynomial_sum(coefficients: list | np.ndarray, t: np.ndarray):
     """
@@ -179,7 +187,7 @@ def mean_longitudes(
     """
     if MEEUS:
         # convert from MJD to days relative to 2000-01-01T12:00:00
-        T = MJD - 51544.5
+        T = MJD - _mjd_j2000
         # mean longitude of moon
         lunar_longitude = np.array([218.3164591, 13.17639647754579,
             -9.9454632e-13, 3.8086292e-20, -8.6184958e-27])
@@ -197,10 +205,10 @@ def mean_longitudes(
             1.55628359e-12, 4.390675353e-20, -9.26940435e-27])
         N = polynomial_sum(lunar_node, T)
         # mean longitude of solar perigee (Simon et al., 1994)
-        PP = 282.94 + 1.7192 * T / 36525.0
+        PP = 282.94 + (1.7192 * T)/_century
     elif ASTRO5:
         # convert from MJD to centuries relative to 2000-01-01T12:00:00
-        T = (MJD - 51544.5)/36525.0
+        T = (MJD - _mjd_j2000)/_century
         # mean longitude of moon (p. 338)
         lunar_longitude = np.array([218.3164477, 481267.88123421, -1.5786e-3,
              1.855835e-6, -1.53388e-8])
@@ -299,7 +307,7 @@ def doodson_arguments(
     # degrees to radians
     dtr = np.pi/180.0
     # convert from MJD to centuries relative to 2000-01-01T12:00:00
-    T = (MJD - 51544.5)/36525.0
+    T = (MJD - _mjd_j2000)/_century
     # hour of the day
     hour = np.mod(MJD, 1)*24.0
     # calculate Doodson phase angles
@@ -383,7 +391,7 @@ def delaunay_arguments(MJD: np.ndarray):
     # arcseconds to radians
     atr = np.pi/648000.0
     # convert from MJD to centuries relative to 2000-01-01T12:00:00
-    T = (MJD - 51544.5)/36525.0
+    T = (MJD - _mjd_j2000)/_century
     # 360 degrees
     circle = 1296000
     # mean anomaly of the moon (arcseconds)
@@ -436,7 +444,7 @@ def mean_obliquity(MJD: np.ndarray):
     # arcseconds to radians
     atr = np.pi/648000.0
     # convert from MJD to centuries relative to 2000-01-01T12:00:00
-    T = (MJD - 51544.5)/36525.0
+    T = (MJD - _mjd_j2000)/_century
     # mean obliquity of the ecliptic (arcseconds)
     epsilon0 = np.array([84381.406, -46.836769, -1.831e-4,
         2.00340e-4, -5.76e-07, -4.34e-08])
@@ -794,9 +802,9 @@ def gast(T: float | np.ndarray):
         <https://iers-conventions.obspm.fr/content/tn36.pdf>`_
     """
     # create timescale from centuries relative to 2000-01-01T12:00:00
-    ts = timescale.time.Timescale(MJD=T*36525.0 + 51544.5)
+    ts = timescale.time.Timescale(MJD=T*_century + _mjd_j2000)
     # convert dynamical time to modified Julian days
-    MJD = ts.tt - 2400000.5
+    MJD = ts.tt - _jd_mjd
     # estimate the mean obliquity
     epsilon = mean_obliquity(MJD)
     # estimate the nutation in longitude and obliquity
@@ -838,9 +846,9 @@ def itrs(T: float | np.ndarray):
         <https://iers-conventions.obspm.fr/content/tn36.pdf>`_
     """
     # create timescale from centuries relative to 2000-01-01T12:00:00
-    ts = timescale.time.Timescale(MJD=T*36525.0 + 51544.5)
+    ts = timescale.time.Timescale(MJD=T*_century + _mjd_j2000)
     # convert dynamical time to modified Julian days
-    MJD = ts.tt - 2400000.5
+    MJD = ts.tt - _jd_mjd
     # estimate the mean obliquity
     epsilon = mean_obliquity(MJD)
     # estimate the nutation in longitude and obliquity
@@ -892,7 +900,7 @@ def _eqeq_complement(T: float | np.ndarray):
         <https://iers-conventions.obspm.fr/content/tn36.pdf>`_
     """
     # create timescale from centuries relative to 2000-01-01T12:00:00
-    ts = timescale.time.Timescale(MJD=T*36525.0 + 51544.5)
+    ts = timescale.time.Timescale(MJD=T*_century + _mjd_j2000)
     # get the fundamental arguments in radians
     fa = np.zeros((14, len(ts)))
     # mean anomaly of the moon (arcseconds)
@@ -988,9 +996,9 @@ def _nutation_angles(T: float | np.ndarray):
         <https://iers-conventions.obspm.fr/content/tn36.pdf>`_
     """
     # create timescale from centuries relative to 2000-01-01T12:00:00
-    ts = timescale.time.Timescale(MJD=T*36525.0 + 51544.5)
+    ts = timescale.time.Timescale(MJD=T*_century + _mjd_j2000)
     # convert dynamical time to modified Julian days
-    MJD = ts.tt - 2400000.5
+    MJD = ts.tt - _jd_mjd
     # get the fundamental arguments in radians
     l, lp, F, D, Om = delaunay_arguments(MJD)
     # non-polynomial terms in the equation of the equinoxes
@@ -1064,7 +1072,7 @@ def _polar_motion_matrix(T: float | np.ndarray):
     # arcseconds to radians
     atr = np.pi/648000.0
     # convert to MJD from centuries relative to 2000-01-01T12:00:00
-    MJD = T*36525.0 + 51544.5
+    MJD = T*_century + _mjd_j2000
     sprime = -4.7e-5*T
     px, py = timescale.eop.iers_polar_motion(MJD)
     # calculate the rotation matrices

diff --git a/pyTMD/compute.py b/pyTMD/compute.py
@@ -61,6 +61,7 @@
 
 UPDATE HISTORY:
     Updated 07/2024: assert that data type is a known value
+        make number of days to convert JD to MJD a variable
     Updated 06/2024: use np.clongdouble instead of np.longcomplex
     Updated 04/2024: use wrapper to importlib for optional dependencies
     Updated 02/2024: changed class name for ellipsoid parameters to datum
@@ -125,6 +126,9 @@
 # attempt imports
 pyproj = pyTMD.utilities.import_dependency('pyproj')
 
+# number of days between the Julian day epoch and MJD
+_jd_mjd = 2400000.5
+
 # PURPOSE: wrapper function for computing values
 def corrections(
         x: np.ndarray, y: np.ndarray, delta_time: np.ndarray,
@@ -828,7 +832,7 @@ def LPT_displacements(
             epoch=EPOCH, standard=TIME)
 
     # convert dynamic time to Modified Julian Days (MJD)
-    MJD = ts.tt - 2400000.5
+    MJD = ts.tt - _jd_mjd
     # convert Julian days to calendar dates
     Y,M,D,h,m,s = timescale.time.convert_julian(ts.tt, format='tuple')
     # calculate time in year-decimal format
@@ -1008,7 +1012,7 @@ def OPT_displacements(
             epoch=EPOCH, standard=TIME)
 
     # convert dynamic time to Modified Julian Days (MJD)
-    MJD = ts.tt - 2400000.5
+    MJD = ts.tt - _jd_mjd
     # convert Julian days to calendar dates
     Y,M,D,h,m,s = timescale.time.convert_julian(ts.tt, format='tuple')
     # calculate time in year-decimal format

diff --git a/pyTMD/predict.py b/pyTMD/predict.py
@@ -21,6 +21,7 @@
 
 UPDATE HISTORY:
     Updated 07/2024: use normalize_angle from pyTMD astro module
+        make number of days to convert tide time to MJD a variable
     Updated 02/2024: changed class name for ellipsoid parameters to datum
     Updated 01/2024: moved minor arguments calculation into new function
         moved constituent parameters function from predict to arguments
@@ -52,6 +53,9 @@
 import pyTMD.astro
 from pyTMD.crs import datum
 
+# number of days between MJD and the tide epoch (1992-01-01T00:00:00)
+_mjd_tide = 48622.0
+
 # PURPOSE: Predict tides at single times
 def map(t: float | np.ndarray,
         hc: np.ndarray,
@@ -93,7 +97,7 @@ def map(t: float | np.ndarray,
     npts, nc = np.shape(hc)
     # load the nodal corrections
     # convert time to Modified Julian Days (MJD)
-    pu, pf, G = pyTMD.arguments.arguments(t + 48622.0,
+    pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide,
         constituents,
         deltat=deltat,
         corrections=corrections
@@ -159,7 +163,7 @@ def drift(t: float | np.ndarray,
     nt = len(t)
     # load the nodal corrections
     # convert time to Modified Julian Days (MJD)
-    pu, pf, G = pyTMD.arguments.arguments(t + 48622.0,
+    pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide,
         constituents,
         deltat=deltat,
         corrections=corrections
@@ -225,7 +229,7 @@ def time_series(t: float | np.ndarray,
     nt = len(t)
     # load the nodal corrections
     # convert time to Modified Julian Days (MJD)
-    pu, pf, G = pyTMD.arguments.arguments(t + 48622.0,
+    pu, pf, G = pyTMD.arguments.arguments(t + _mjd_tide,
         constituents,
         deltat=deltat,
         corrections=corrections
@@ -369,7 +373,7 @@ def infer_minor(
 
     # load the nodal corrections for minor constituents
     # convert time to Modified Julian Days (MJD)
-    pu, pf, G = pyTMD.arguments.minor_arguments(t + 48622.0,
+    pu, pf, G = pyTMD.arguments.minor_arguments(t + _mjd_tide,
         deltat=kwargs['deltat'],
         corrections=kwargs['corrections']
     )
@@ -535,7 +539,7 @@ def solid_earth_tide(
     # number of input coordinates
     nt = len(np.atleast_1d(t))
     # convert time to Modified Julian Days (MJD)
-    MJD = t + 48622.0
+    MJD = t + _mjd_tide
     # scalar product of input coordinates with sun/moon vectors
     radius = np.sqrt(XYZ[:,0]**2 + XYZ[:,1]**2 + XYZ[:,2]**2)
     solar_radius = np.sqrt(SXYZ[:,0]**2 + SXYZ[:,1]**2 + SXYZ[:,2]**2)

diff --git a/test/test_solid_earth.py b/test/test_solid_earth.py
@@ -96,6 +96,8 @@ def test_frequency_dependence_diurnal():
     T = (MJD - 51544.5)/36525.0
     T_expected = 0.1059411362080767
     assert np.isclose(T_expected, T)
+    T_test = (MJD - pyTMD.astro._mjd_j2000)/pyTMD.astro._century
+    assert np.isclose(T_expected, T_test)
     # expected results
     dx_expected = 0.4193085327321284701e-2
     dy_expected = 0.1456681241014607395e-2
@@ -116,6 +118,8 @@ def test_frequency_dependence_long_period():
     T = (MJD - 51544.5)/36525.0
     T_expected = 0.1059411362080767
     assert np.isclose(T_expected, T)
+    T_test = (MJD - pyTMD.astro._mjd_j2000)/pyTMD.astro._century
+    assert np.isclose(T_expected, T_test)
     # expected results
     dx_expected = -0.9780962849562107762e-4
     dy_expected = -0.2236349699932734273e-4
@@ -148,6 +152,8 @@ def test_fundamental_arguments():
     # convert to MJD from centuries relative to 2000-01-01T12:00:00
     MJD = T*36525.0 + 51544.5
     assert np.isclose(MJD, 54465)
+    MJD_test = T*pyTMD.astro._century + pyTMD.astro._mjd_j2000
+    assert np.isclose(MJD, MJD_test)
     L_expected = 2.291187512612069099
     LP_expected = 6.212931111003726414
     F_expected = 3.658025792050572989