feat: add properties for Doodson and Cartwright numbers to `constitue…

…nts` class

pyTMD/arguments.py

@@ -502,6 +502,7 @@ def arguments(
     G = np.zeros((nt,nc))
     for i,c in enumerate(constituents):
         # map between given constituents and supported in tidal program
+        assert c.lower() in cindex, f'Unsupported constituent {c.lower()}'
         j, = [j for j,val in enumerate(cindex) if (val == c.lower())]
         pu[:,i] = u[:,j]*dtr
         pf[:,i] = f[:,j]
@@ -707,6 +708,8 @@ def doodson_number(
     # set default keyword arguments
     kwargs.setdefault('corrections', 'OTIS')
     kwargs.setdefault('formalism', 'Doodson')
+    # validate inputs
+    assert kwargs['formalism'] in ('Cartwright', 'Doodson'), 'Unknown formalism'
 
     # constituents array (not all are included in tidal program)
     cindex = ['sa', 'ssa', 'mm', 'msf', 'mf', 'mt', 'alpha1', '2q1', 'sigma1',
@@ -719,6 +722,7 @@ def doodson_number(
     coefficients = _arguments_table(**kwargs)
     if isinstance(constituents, str):
         # map between given constituents and supported in tidal program
+        assert constituents.lower() in cindex, f'Unsupported constituent'
         j, = [j for j,val in enumerate(cindex) if (val == constituents.lower())]
         # extract identifier in formalism
         if (kwargs['formalism'] == 'Cartwright'):
@@ -733,6 +737,7 @@ def doodson_number(
         # for each input constituent
         for i,c in enumerate(constituents):
             # map between given constituents and supported in tidal program
+            assert c.lower() in cindex, f'Unsupported constituent {c.lower()}'
             j, = [j for j,val in enumerate(cindex) if (val == c.lower())]
             # convert from coefficients to Doodson number
             if (kwargs['formalism'] == 'Cartwright'):

pyTMD/io/constituents.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 constituents.py
-Written by Tyler Sutterley (08/2023)
+Written by Tyler Sutterley (01/2024)
 Basic tide model constituent class
 
 PYTHON DEPENDENCIES:
@@ -10,6 +10,7 @@
         https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
 
 UPDATE HISTORY:
+    Updated 01/2024: added properties for Doodson and Cartwright numbers
     Updated 08/2023: added default for printing constituent class
     Updated 07/2023: output constituent from get and pop as copy
     Updated 03/2023: add basic variable typing to function inputs
@@ -19,6 +20,7 @@
 
 import copy
 import numpy as np
+import pyTMD.arguments
 
 class constituents:
     """
@@ -152,6 +154,19 @@ def phase(self, field: str):
         ph.data[ph.mask] = ph.fill_value
         return ph
 
+    @property
+    def doodson_number(self):
+        """Doodson number of constituents
+        """
+        return [pyTMD.arguments.doodson_number(f) for f in self.fields]
+
+    @property
+    def cartwright_number(self):
+        """Cartwright numbers of constituents
+        """
+        return [pyTMD.arguments.doodson_number(f, formalism='Cartwright')
+                for f in self.fields]
+
     def __str__(self):
         """String representation of the ``constituents`` object
         """

test/test_atlas_read.py

@@ -38,6 +38,7 @@
 import pyTMD.io
 import pyTMD.time
 import pyTMD.io.model
+import pyTMD.arguments
 import pyTMD.utilities
 
 # current file path
@@ -197,6 +198,13 @@ def test_compare_TPXO9_v2(METHOD):
     # calculate complex form of constituent oscillation
     hc2 = amp2*np.exp(-1j*ph2*np.pi/180.0)
 
+    # expected Doodson numbers for constituents
+    exp = {}
+    exp['m2'] = 255.555
+    exp['s2'] = 273.555
+    exp['o1'] = 145.555
+    exp['k1'] = 165.555
+
     # will verify differences between model outputs are within tolerance
     eps = np.finfo(np.float16).eps
     # calculate differences between methods
@@ -206,6 +214,11 @@ def test_compare_TPXO9_v2(METHOD):
         # calculate difference in amplitude and phase
         difference = hc1[:,i] - hc2[:,i]
         assert np.all(np.abs(difference) <= eps)
+        # verify doodson numbers
+        assert (constituents.doodson_number[i] == exp[cons])
+        # verify cartwright numbers
+        assert np.all(constituents.cartwright_number[i] ==
+                pyTMD.arguments._from_doodson_number(exp[cons]))
 
 # parameterize interpolation method
 @pytest.mark.parametrize("METHOD", ['bilinear'])

test/test_fes_predict.py

@@ -192,13 +192,36 @@ def test_compare_FES2014(METHOD):
     # calculate complex form of constituent oscillation
     hc2 = amp2*np.exp(-1j*ph2*np.pi/180.0)
 
+    # expected Doodson numbers for constituents
+    exp = {}
+    exp['2n2'] = 235.755
+    exp['k1'] = 165.555
+    exp['k2'] = 275.555
+    exp['m2'] = 255.555
+    exp['m4'] = 455.555
+    exp['mf'] = 75.555
+    exp['mm'] = 65.455
+    exp['msqm'] = 93.555
+    exp['mtm'] = 85.455
+    exp['n2'] = 245.655
+    exp['o1'] = 145.555
+    exp['p1'] = 163.555
+    exp['q1'] = 135.655
+    exp['s1'] = 164.555
+    exp['s2'] = 273.555
+
     # will verify differences between model outputs are within tolerance
     eps = np.finfo(np.float16).eps
     # calculate differences between methods
     for i, cons in enumerate(c):
         # calculate difference in amplitude and phase
         difference = hc1[:,i] - hc2[:,i]
         assert np.all(np.abs(difference) <= eps)
+        # verify doodson numbers
+        assert (constituents.doodson_number[i] == exp[cons])
+        # verify cartwright numbers
+        assert np.all(constituents.cartwright_number[i] ==
+                pyTMD.arguments._from_doodson_number(exp[cons]))
 
     # validate iteration within constituents class
     for field, hc in constituents: