diff --git a/pytides2/constituent.py b/pytides2/constituent.py index 6c511d8..878a9f2 100644 --- a/pytides2/constituent.py +++ b/pytides2/constituent.py @@ -140,6 +140,7 @@ def f(self, a: Dict[str, AstronomicalParameter]) -> Union[Any, float, int]: _Z0 = BaseConstituent(name="Z0", xdo="Z ZZZ ZZZ", u=nc.u_zero, f=nc.f_unity) _Sa = BaseConstituent(name="Sa", xdo="Z ZAZ ZZZ", u=nc.u_zero, f=nc.f_unity) _Ssa = BaseConstituent(name="Ssa", xdo="Z ZBZ ZZZ", u=nc.u_zero, f=nc.f_unity) +_MSm = BaseConstituent(name="MSm", xdo="Z AXA ZZZ", u=nc.u_zero, f=nc.f_Mm) _Mm = BaseConstituent(name="Mm", xdo="Z AZY ZZZ", u=nc.u_zero, f=nc.f_Mm) _Mf = BaseConstituent(name="Mf", xdo="Z BZZ ZZZ", u=nc.u_Mf, f=nc.f_Mf) @@ -159,6 +160,22 @@ def f(self, a: Dict[str, AstronomicalParameter]) -> Union[Any, float, int]: _S1 = BaseConstituent(name="S1", xdo="A AYZ ZZZ", u=nc.u_zero, f=nc.f_unity) _OO1 = BaseConstituent(name="OO1", xdo="A CZZ ZZY", u=nc.u_OO1, f=nc.f_OO1) +# Additional minor diurnals. The lunar members are grouped with the major +# constituent sharing their node-factor structure (O1 for the lower-frequency +# group, J1 for the upper group, OO1 for the highest); the solar members +# (pi1, psi1) require no nodal correction. +_alpha1 = BaseConstituent(name="alpha1", xdo="A VBA ZZA", u=nc.u_O1, f=nc.f_O1) +_sigma1 = BaseConstituent(name="sigma1", xdo="A WBZ ZZA", u=nc.u_O1, f=nc.f_O1) +_tau1 = BaseConstituent(name="tau1", xdo="A YBZ ZZA", u=nc.u_O1, f=nc.f_O1) +_beta1 = BaseConstituent(name="beta1", xdo="A ZXA ZZA", u=nc.u_O1, f=nc.f_O1) +_NO1 = BaseConstituent(name="NO1", xdo="A ZZA ZZA", u=nc.u_J1, f=nc.f_J1) +_chi1 = BaseConstituent(name="chi1", xdo="A ZBY ZZA", u=nc.u_J1, f=nc.f_J1) +_pi1 = BaseConstituent(name="pi1", xdo="A AWZ ZAA", u=nc.u_zero, f=nc.f_unity) +_psi1 = BaseConstituent(name="psi1", xdo="A AAZ ZYY", u=nc.u_zero, f=nc.f_unity) +_phi1 = BaseConstituent(name="phi1", xdo="A ABZ ZZY", u=nc.u_J1, f=nc.f_J1) +_theta1 = BaseConstituent(name="theta1", xdo="A BXA ZZA", u=nc.u_J1, f=nc.f_J1) +_upsilon1 = BaseConstituent(name="upsilon1", xdo="A DZY ZZY", u=nc.u_OO1, f=nc.f_OO1) + # Semi-Diurnals _2N2 = BaseConstituent(name="2N2", xdo="B XZB ZZZ", u=nc.u_M2, f=nc.f_M2) _N2 = BaseConstituent(name="N2", xdo="B YZA ZZZ", u=nc.u_M2, f=nc.f_M2) @@ -171,6 +188,12 @@ def f(self, a: Dict[str, AstronomicalParameter]) -> Union[Any, float, int]: _R2 = BaseConstituent(name="R2", xdo="B BYZ ZYB", u=nc.u_zero, f=nc.f_unity) _K2 = BaseConstituent(name="K2", xdo="B BZZ ZZZ", u=nc.u_K2, f=nc.f_K2) +# Additional minor semi-diurnals. OQ2 is a lunar elliptic line (M2 node +# factor); H1 and H2 are solar lines requiring no nodal correction. +_OQ2 = BaseConstituent(name="OQ2", xdo="B WZC ZZZ", u=nc.u_M2, f=nc.f_M2) +_H1 = BaseConstituent(name="H1", xdo="B ZYZ ZAZ", u=nc.u_zero, f=nc.f_unity) +_H2 = BaseConstituent(name="H2", xdo="B ZAZ ZYZ", u=nc.u_zero, f=nc.f_unity) + # Third-Diurnals _M3 = BaseConstituent( name="M3", xdo="C ZZZ ZZZ", u=lambda a: nc.u_Modd(a, 3), f=lambda a: nc.f_Modd(a, 3) @@ -183,30 +206,55 @@ def f(self, a: Dict[str, AstronomicalParameter]) -> Union[Any, float, int]: # Diurnal _2Q1 = CompoundConstituent(name="2Q1", members=[(_N2, 1), (_J1, -1)]) _rho1 = CompoundConstituent(name="rho1", members=[(_nu2, 1), (_K1, -1)]) +_SO1 = CompoundConstituent(name="SO1", members=[(_S2, 1), (_O1, -1)]) # Semi-Diurnal _mu2 = CompoundConstituent(name="mu2", members=[(_M2, 2), (_S2, -1)]) # 2MS2 _2SM2 = CompoundConstituent(name="2SM2", members=[(_S2, 2), (_M2, -1)]) +_epsilon2 = CompoundConstituent( # MNS2 + name="epsilon2", members=[(_M2, 1), (_N2, 1), (_S2, -1)] +) +_MKS2 = CompoundConstituent(name="MKS2", members=[(_M2, 1), (_K2, 1), (_S2, -1)]) +_MSN2 = CompoundConstituent(name="MSN2", members=[(_M2, 1), (_S2, 1), (_N2, -1)]) +_eta2 = CompoundConstituent(name="eta2", members=[(_K1, 1), (_J1, 1)]) # KJ2 # Third-Diurnal _2MK3 = CompoundConstituent(name="2MK3", members=[(_M2, 1), (_O1, 1)]) _MK3 = CompoundConstituent(name="MK3", members=[(_M2, 1), (_K1, 1)]) +_MO3 = CompoundConstituent(name="MO3", members=[(_M2, 1), (_O1, 1)]) +_SO3 = CompoundConstituent(name="SO3", members=[(_S2, 1), (_O1, 1)]) +_SK3 = CompoundConstituent(name="SK3", members=[(_S2, 1), (_K1, 1)]) # Quarter-Diurnal _MN4 = CompoundConstituent(name="MN4", members=[(_M2, 1), (_N2, 1)]) _M4 = CompoundConstituent(name="M4", members=[(_M2, 2)]) _MS4 = CompoundConstituent(name="MS4", members=[(_M2, 1), (_S2, 1)]) _S4 = CompoundConstituent(name="S4", members=[(_S2, 2)]) +_SN4 = CompoundConstituent(name="SN4", members=[(_S2, 1), (_N2, 1)]) +_MK4 = CompoundConstituent(name="MK4", members=[(_M2, 1), (_K2, 1)]) +_SK4 = CompoundConstituent(name="SK4", members=[(_S2, 1), (_K2, 1)]) + +# Fifth-Diurnal +_2MK5 = CompoundConstituent(name="2MK5", members=[(_M2, 2), (_K1, 1)]) +_2SK5 = CompoundConstituent(name="2SK5", members=[(_S2, 2), (_K1, 1)]) # Sixth-Diurnal _M6 = CompoundConstituent(name="M6", members=[(_M2, 3)]) _S6 = CompoundConstituent(name="S6", members=[(_S2, 3)]) +_2MN6 = CompoundConstituent(name="2MN6", members=[(_M2, 2), (_N2, 1)]) +_2MS6 = CompoundConstituent(name="2MS6", members=[(_M2, 2), (_S2, 1)]) +_2MK6 = CompoundConstituent(name="2MK6", members=[(_M2, 2), (_K2, 1)]) +_2SM6 = CompoundConstituent(name="2SM6", members=[(_M2, 1), (_S2, 2)]) +_MSK6 = CompoundConstituent(name="MSK6", members=[(_M2, 1), (_S2, 1), (_K2, 1)]) + +# Seventh-Diurnal +_3MK7 = CompoundConstituent(name="3MK7", members=[(_M2, 3), (_K1, 1)]) # Eighth-Diurnals _M8 = CompoundConstituent(name="M8", members=[(_M2, 4)]) - +# NOAA set (37 tidal constituents) noaa = [ _M2, _S2, @@ -246,3 +294,77 @@ def f(self, a: Dict[str, AstronomicalParameter]) -> Union[Any, float, int]: _M8, _MS4, ] + +# Extended 67-constituent set (the NOAA set above augmented with the additional +# minor and shallow-water constituents), ordered by increasing frequency. +# list provided by Pengcheng Wang, obtained by setting Rayleigh criteria to 0.8 +# note M1 in NOAA set is replaced by NO1 here. +extended = [ + _Sa, + _Ssa, + _MSm, + _Mm, + _MSF, + _Mf, + _alpha1, + _2Q1, + _sigma1, + _Q1, + _rho1, + _O1, + _tau1, + _beta1, + _NO1, + _chi1, + _pi1, + _P1, + _S1, + _K1, + _psi1, + _phi1, + _theta1, + _J1, + _SO1, + _OO1, + _upsilon1, + _OQ2, + _epsilon2, + _2N2, + _mu2, + _N2, + _nu2, + _H1, + _M2, + _H2, + _MKS2, + _lambda2, + _L2, + _T2, + _S2, + _R2, + _K2, + _MSN2, + _eta2, + _MO3, + _M3, + _SO3, + _MK3, + _SK3, + _MN4, + _M4, + _SN4, + _MS4, + _MK4, + _S4, + _SK4, + _2MK5, + _2SK5, + _2MN6, + _M6, + _2MS6, + _2MK6, + _2SM6, + _MSK6, + _3MK7, + _M8, +] diff --git a/pytides2/nodal_corrections.py b/pytides2/nodal_corrections.py index 728c748..876e3ae 100644 --- a/pytides2/nodal_corrections.py +++ b/pytides2/nodal_corrections.py @@ -103,7 +103,7 @@ def f_K2(a: Dict[str, AstronomicalParameter]) -> Union[float, Any]: sinsqIcos2nu_mean: float = np.sin(omega) ** 2 * (1 - 3 / 2.0 * np.sin(i) ** 2) mean = 0.5023 * sinsqIcos2nu_mean + 0.0365 return ( - 0.2533 * np.sin(I) ** 4 + 0.0367 * np.sin(I) ** 2 * np.cos(2 * nu) + 0.0013 + 0.2523 * np.sin(I) ** 4 + 0.0367 * np.sin(I) ** 2 * np.cos(2 * nu) + 0.0013 ) ** (0.5) / mean diff --git a/test/data/Tidal_constituents_67.csv b/test/data/Tidal_constituents_67.csv new file mode 100644 index 0000000..c27976f --- /dev/null +++ b/test/data/Tidal_constituents_67.csv @@ -0,0 +1,68 @@ +Name,Frequency +SA ,0.000114074 +SSA ,0.000228159 +MSM ,0.001309781 +MM ,0.001512152 +MSF ,0.002821933 +MF ,0.003050092 +ALP1,0.03439657 +2Q1 ,0.035706351 +SIG1,0.035908722 +Q1 ,0.037218502 +RHO1,0.037420874 +O1 ,0.038730654 +TAU1,0.038958814 +BET1,0.040040435 +NO1 ,0.040268594 +CHI1,0.040470966 +PI1 ,0.041438513 +P1 ,0.041552587 +S1 ,0.041666672 +K1 ,0.041780746 +PSI1,0.04189482 +PHI1,0.042008905 +THE1,0.043090527 +J1 ,0.043292898 +SO1 ,0.044602679 +OO1 ,0.044830838 +UPS1,0.04634299 +OQ2 ,0.075974945 +EPS2,0.076177316 +2N2 ,0.077487097 +MU2 ,0.077689468 +N2 ,0.078999249 +NU2 ,0.07920162 +H1 ,0.080397327 +M2 ,0.080511401 +H2 ,0.080625475 +MKS2,0.08073956 +LDA2,0.081821181 +L2 ,0.082023553 +T2 ,0.083219259 +S2 ,0.083333333 +R2 ,0.083447407 +K2 ,0.083561492 +MSN2,0.084845485 +ETA2,0.085073644 +MO3 ,0.119242055 +M3 ,0.120767101 +SO3 ,0.122063988 +MK3 ,0.122292147 +SK3 ,0.12511408 +MN4 ,0.159510649 +M4 ,0.161022801 +SN4 ,0.162332582 +MS4 ,0.163844734 +MK4 ,0.164072893 +S4 ,0.166666667 +SK4 ,0.166894826 +2MK5,0.202803548 +2SK5,0.208447413 +2MN6,0.24002205 +M6 ,0.241534202 +2MS6,0.244356135 +2MK6,0.244584294 +2SM6,0.247178067 +MSK6,0.247406226 +3MK7,0.283314948 +M8 ,0.322045603 \ No newline at end of file diff --git a/test/test_pytides2.py b/test/test_pytides2.py index 2b75c99..807e464 100644 --- a/test/test_pytides2.py +++ b/test/test_pytides2.py @@ -1,4 +1,5 @@ -from datetime import datetime +import csv +from datetime import datetime, timedelta import numpy as np @@ -7,6 +8,93 @@ def test_import() -> None: from pytides2 import astro, constituent, nodal_corrections, tide +def test_extended_constituents() -> None: + from pytides2 import constituent + from pytides2.astro import astro + + extended = constituent.extended + # The extended set reproduces the bundled 67-constituent table. + assert len(extended) == 67 + assert len({c.name for c in extended}) == 67 + + # pytides spells out Greek letters in lower case where the table uses an + # all-upper-case abbreviation. + alias = { + "ALP1": "alpha1", + "SIG1": "sigma1", + "BET1": "beta1", + "THE1": "theta1", + "UPS1": "upsilon1", + "EPS2": "epsilon2", + "LDA2": "lambda2", + "MSM": "MSm", + } + by_name = {c.name.upper(): c for c in extended} + a = astro(datetime(2010, 6, 15, 7)) + with open("test/data/Tidal_constituents_67.csv", encoding="utf-8-sig") as fh: + reader = csv.reader(fh) + next(reader) + rows = [(name.strip(), float(freq)) for name, freq in reader] + + assert len(rows) == 67 + for name, frequency in rows: + c = by_name[alias.get(name, name).upper()] + # speed() is in degrees/hour; the table lists cycles/hour. + assert abs(c.speed(a) / 360.0 - frequency) < 1e-6, name + + +def test_Tide_noaa() -> None: + from pytides2 import constituent + from pytides2.tide import Tide + + by_name = {c.name: c for c in constituent.noaa} + amplitudes = {"M2": 1.5, "S2": 0.5, "K1": 0.4, "O1": 0.3, "N2": 0.2} + phases = {"M2": 45.0, "S2": 90.0, "K1": 120.0, "O1": 200.0, "N2": 300.0} + + synthetic = Tide( + constituents=[by_name[n] for n in amplitudes], + amplitudes=[amplitudes[n] for n in amplitudes], + phases=[phases[n] for n in amplitudes], + ) + t0 = datetime(2015, 1, 1) + times = [t0 + timedelta(hours=h) for h in range(24 * 120)] + heights = synthetic.at(times) + + # Fitting with the NOAA set should recover the input constituents. + tide, _ = Tide.decompose(heights, np.array(times), constituents=constituent.noaa) + fitted = {m["constituent"].name: m for m in tide.model} + for name in amplitudes: + assert abs(fitted[name]["amplitude"] - amplitudes[name]) < 1e-3 + assert abs(fitted[name]["phase"] - phases[name]) < 0.5 + + +def test_Tide_extended() -> None: + from pytides2 import constituent + from pytides2.tide import Tide + + by_name = {c.name: c for c in constituent.extended} + amplitudes = {"M2": 1.5, "S2": 0.5, "K1": 0.4, "O1": 0.3, "N2": 0.2} + phases = {"M2": 45.0, "S2": 90.0, "K1": 120.0, "O1": 200.0, "N2": 300.0} + + synthetic = Tide( + constituents=[by_name[n] for n in amplitudes], + amplitudes=[amplitudes[n] for n in amplitudes], + phases=[phases[n] for n in amplitudes], + ) + t0 = datetime(2015, 1, 1) + times = [t0 + timedelta(hours=h) for h in range(24 * 120)] + heights = synthetic.at(times) + + # Fitting with the full extended set should recover the input constituents. + tide, _ = Tide.decompose( + heights, np.array(times), constituents=constituent.extended + ) + fitted = {m["constituent"].name: m for m in tide.model} + for name in amplitudes: + assert abs(fitted[name]["amplitude"] - amplitudes[name]) < 1e-3 + assert abs(fitted[name]["phase"] - phases[name]) < 0.5 + + def test_Tide() -> None: from pytides2.tide import Tide @@ -29,9 +117,11 @@ def test_Tide() -> None: assert t.model[0]["amplitude"] == 178.94736842105263 assert t.model[0][2] == 0.0 assert len(t.model) == 33 - assert ( - t.formzahl == 0.3691950696609845 - or t.formzahl == 0.38101230038937306 - or t.formzahl == 0.38376753628038296 + # The exact value drifts slightly with the numpy/scipy/BLAS version (the + # least-squares solver does not fully converge on this short record), so + # compare against the known results with a tolerance rather than exactly. + assert any( + abs(t.formzahl - v) < 1e-3 + for v in (0.3691950696609845, 0.38101230038937306, 0.38376753628038296) ) assert t.type == "mixed (semidiurnal)"