From a57d4daae781dc20fdb987f3633a28764f42b4c0 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Sun, 1 Dec 2024 11:08:12 -0800 Subject: [PATCH] refactor: move body tide Love/Shida numbers to `arguments` --- doc/source/api_reference/arguments.rst | 2 + doc/source/api_reference/predict.rst | 2 - pyTMD/arguments.py | 101 ++++++++++++++++++++++++ pyTMD/predict.py | 105 +------------------------ test/test_love_numbers.py | 6 +- 5 files changed, 109 insertions(+), 107 deletions(-) diff --git a/doc/source/api_reference/arguments.rst b/doc/source/api_reference/arguments.rst index 0f238502..c9519c7f 100644 --- a/doc/source/api_reference/arguments.rst +++ b/doc/source/api_reference/arguments.rst @@ -36,6 +36,8 @@ Calling Sequence .. autofunction:: pyTMD.arguments._constituent_parameters +.. autofunction:: pyTMD.arguments._love_numbers + .. autofunction:: pyTMD.arguments._parse_tide_potential_table .. autofunction:: pyTMD.arguments._to_doodson_number diff --git a/doc/source/api_reference/predict.rst b/doc/source/api_reference/predict.rst index 484091d1..43d1e67a 100644 --- a/doc/source/api_reference/predict.rst +++ b/doc/source/api_reference/predict.rst @@ -40,8 +40,6 @@ Calling Sequence .. autofunction:: pyTMD.predict._infer_long_period -.. autofunction:: pyTMD.predict._body_tide_love_numbers - .. autofunction:: pyTMD.predict.equilibrium_tide .. autofunction:: pyTMD.predict.load_pole_tide diff --git a/pyTMD/arguments.py b/pyTMD/arguments.py index 5202c5a2..9937ef50 100755 --- a/pyTMD/arguments.py +++ b/pyTMD/arguments.py @@ -40,6 +40,7 @@ UPDATE HISTORY: Updated 11/2024: allow variable case for Doodson number formalisms fix species in constituent parameters for complex tides + move body tide Love/Shida numbers from predict module Updated 10/2024: can convert Doodson numbers formatted as strings update Doodson number conversions to follow Cartwright X=10 convention add function to parse Cartwright/Tayler/Edden tables @@ -93,6 +94,7 @@ "_arguments_table", "_minor_table", "_constituent_parameters", + "_love_numbers", "_parse_tide_potential_table", "_to_doodson_number", "_to_extended_doodson", @@ -1462,6 +1464,105 @@ def _constituent_parameters(c: str, **kwargs): # return the values for the constituent return (amplitude, phase, omega, alpha, species) +def _love_numbers( + omega: np.ndarray, + model: str = 'PREM' + ): + """ + Compute the body tide Love/Shida numbers for a given + frequency [1]_ [2]_ [3]_ + + Parameters + ---------- + omega: np.ndarray + angular frequency (radians per second) + model: str, default 'PREM' + Earth model to use for Love numbers + + - '1066A' + - 'PEM-C' + - 'C2' + - 'PREM' + + Returns + ------- + h2: float + Degree-2 Love number of vertical displacement + k2: float + Degree-2 Love number of gravitational potential + l2: float + Degree-2 Love (Shida) number of horizontal displacement + + References + ---------- + .. [1] J. M. Wahr, "Body tides on an elliptical, rotating, elastic + and oceanless Earth", *Geophysical Journal of the Royal + Astronomical Society*, 64(3), 677--703, (1981). + `doi: 10.1111/j.1365-246X.1981.tb02690.x + `_ + .. [2] J. M. Wahr and T. Sasao, "A diurnal resonance in the ocean + tide and in the Earth's load response due to the resonant free + `core nutation`", *Geophysical Journal of the Royal Astronomical + Society*, 64(3), 747--765, (1981). + `doi: 10.1111/j.1365-246X.1981.tb02693.x + `_ + .. [3] P. M. Mathews, B. A. Buffett, and I. I. Shapiro, + "Love numbers for diurnal tides: Relation to wobble admittances + and resonance expansions", *Journal of Geophysical Research: + Solid Earth*, 100(B6), 9935--9948, (1995). + `doi: 10.1029/95jb00670 `_ + """ + # free core nutation frequencies (cycles per sidereal day) and + # Love number parameters from Wahr (1981) table 6 + # and Mathews et al. (1995) table 3 + if (model == '1066A'): + fcn = 1.0021714 + h0, h1 = np.array([6.03e-1, -2.46e-3]) + k0, k1 = np.array([2.98e-1, -1.23e-3]) + l0, l1 = np.array([8.42e-2, 7.81e-5]) + elif (model == 'PEM-C'): + fcn = 1.0021771 + h0, h1 = np.array([6.02e-1, -2.46e-3]) + k0, k1 = np.array([2.98e-1, -1.24e-3]) + l0, l1 = np.array([8.39e-2, 7.69e-5]) + elif (model == 'C2'): + fcn = 1.0021844 + h0, h1 = np.array([6.02e-1, -2.45e-3]) + k0, k1 = np.array([2.98e-1, -1.23e-3]) + l0, l1 = np.array([8.46e-2, 7.58e-5]) + elif (model == 'PREM'): + fcn = 1.0023214 + h0, h1 = np.array([5.994e-1, -2.532e-3]) + k0, k1 = np.array([2.962e-1, -1.271e-3]) + l0, l1 = np.array([8.378e-2, 7.932e-5]) + else: + raise ValueError(f'Unknown Earth model: {model}') + # Love numbers for different frequency bands + if (omega > 1e-4): + # tides in the semi-diurnal band + h2 = 0.609 + k2 = 0.302 + l2 = 0.0852 + elif (omega < 2e-5): + # tides in the long period band + h2 = 0.606 + k2 = 0.299 + l2 = 0.0840 + else: + # use resonance formula for tides in the diurnal band + # frequency of the o1 tides (radians/second) + omega_o1, = frequency('o1') + # frequency of free core nutation (radians/second) + omega_fcn = fcn*7292115e-11 + # Love numbers for frequency using equation 4.18 of Wahr (1981) + # (simplification to use only the free core nutation term) + ratio = (omega - omega_o1)/(omega_fcn - omega) + h2 = h0 + h1*ratio + k2 = k0 + k1*ratio + l2 = l0 + l1*ratio + # return the Love numbers for frequency + return (h2, k2, l2) + # Cartwright and Tayler (1971) table with 3rd-degree values _ct1971_table_5 = get_data_path(['data','ct1971_tab5.txt']) # Cartwright and Edden (1973) table with updated values diff --git a/pyTMD/predict.py b/pyTMD/predict.py index 0ffb9026..e466541a 100644 --- a/pyTMD/predict.py +++ b/pyTMD/predict.py @@ -21,6 +21,7 @@ UPDATE HISTORY: Updated 11/2024: use Love numbers for long-period tides when inferring + move body tide Love/Shida numbers to arguments module Updated 10/2024: use PREM as the default Earth model for Love numbers more descriptive error message if cannot infer minor constituents updated calculation of long-period equilibrium tides @@ -80,7 +81,6 @@ "_infer_semi_diurnal", "_infer_diurnal", "_infer_long_period", - "_body_tide_love_numbers", "equilibrium_tide", "load_pole_tide", "ocean_pole_tide", @@ -772,7 +772,7 @@ def _infer_diurnal( if j: j1, = j # Love numbers of degree 2 for constituent - h2, k2, l2 = _body_tide_love_numbers(omajor[i]) + h2, k2, l2 = pyTMD.arguments._love_numbers(omajor[i]) # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # "normalize" tide values @@ -835,7 +835,7 @@ def _infer_diurnal( # sum over the minor tidal constituents of interest for k in minor_indices: # Love numbers of degree 2 for constituent - h2, k2, l2 = _body_tide_love_numbers(omega[k]) + h2, k2, l2 = pyTMD.arguments._love_numbers(omega[k]) # tilt factor: response with respect to the solid earth gamma_2 = (1.0 + k2 - h2) # interpolate from major constituents @@ -1005,105 +1005,6 @@ def _infer_long_period( # return the inferred values return dh -def _body_tide_love_numbers( - omega: np.ndarray, - model: str = 'PREM' - ): - """ - Compute the body tide Love/Shida numbers for a given - frequency [1]_ [2]_ [3]_ - - Parameters - ---------- - omega: np.ndarray - angular frequency (radians per second) - model: str, default 'PREM' - Earth model to use for Love numbers - - - '1066A' - - 'PEM-C' - - 'C2' - - 'PREM' - - Returns - ------- - h2: float - Degree-2 Love number of vertical displacement - k2: float - Degree-2 Love number of gravitational potential - l2: float - Degree-2 Love (Shida) number of horizontal displacement - - References - ---------- - .. [1] J. M. Wahr, "Body tides on an elliptical, rotating, elastic - and oceanless Earth", *Geophysical Journal of the Royal - Astronomical Society*, 64(3), 677--703, (1981). - `doi: 10.1111/j.1365-246X.1981.tb02690.x - `_ - .. [2] J. M. Wahr and T. Sasao, "A diurnal resonance in the ocean - tide and in the Earth's load response due to the resonant free - `core nutation`", *Geophysical Journal of the Royal Astronomical - Society*, 64(3), 747--765, (1981). - `doi: 10.1111/j.1365-246X.1981.tb02693.x - `_ - .. [3] P. M. Mathews, B. A. Buffett, and I. I. Shapiro, - "Love numbers for diurnal tides: Relation to wobble admittances - and resonance expansions", *Journal of Geophysical Research: - Solid Earth*, 100(B6), 9935--9948, (1995). - `doi: 10.1029/95jb00670 `_ - """ - # free core nutation frequencies (cycles per sidereal day) and - # Love number parameters from Wahr (1981) table 6 - # and Mathews et al. (1995) table 3 - if (model == '1066A'): - fcn = 1.0021714 - h0, h1 = np.array([6.03e-1, -2.46e-3]) - k0, k1 = np.array([2.98e-1, -1.23e-3]) - l0, l1 = np.array([8.42e-2, 7.81e-5]) - elif (model == 'PEM-C'): - fcn = 1.0021771 - h0, h1 = np.array([6.02e-1, -2.46e-3]) - k0, k1 = np.array([2.98e-1, -1.24e-3]) - l0, l1 = np.array([8.39e-2, 7.69e-5]) - elif (model == 'C2'): - fcn = 1.0021844 - h0, h1 = np.array([6.02e-1, -2.45e-3]) - k0, k1 = np.array([2.98e-1, -1.23e-3]) - l0, l1 = np.array([8.46e-2, 7.58e-5]) - elif (model == 'PREM'): - fcn = 1.0023214 - h0, h1 = np.array([5.994e-1, -2.532e-3]) - k0, k1 = np.array([2.962e-1, -1.271e-3]) - l0, l1 = np.array([8.378e-2, 7.932e-5]) - else: - raise ValueError(f'Unknown Earth model: {model}') - # Love numbers for different frequency bands - if (omega > 1e-4): - # tides in the semi-diurnal band - h2 = 0.609 - k2 = 0.302 - l2 = 0.0852 - elif (omega < 2e-5): - # tides in the long period band - h2 = 0.606 - k2 = 0.299 - l2 = 0.0840 - else: - # use resonance formula for tides in the diurnal band - # frequency of the o1 tides (radians/second) - omega_o1, = pyTMD.arguments.frequency('o1') - # frequency of free core nutation (radians/second) - omega_fcn = fcn*7292115e-11 - # Love numbers for frequency using equation 4.18 of Wahr (1981) - # (simplification to use only the free core nutation term) - ratio = (omega - omega_o1)/(omega_fcn - omega) - h2 = h0 + h1*ratio - k2 = k0 + k1*ratio - l2 = l0 + l1*ratio - # return the Love numbers for frequency - return (h2, k2, l2) - # PURPOSE: estimate long-period equilibrium tides def equilibrium_tide( t: np.ndarray, diff --git a/test/test_love_numbers.py b/test/test_love_numbers.py index bc393ed3..970a72d1 100644 --- a/test/test_love_numbers.py +++ b/test/test_love_numbers.py @@ -1,15 +1,15 @@ #!/usr/bin/env python u""" -test_love_numbers.py (09/2024) +test_love_numbers.py (11/2024) Verify Love numbers correspond to those from Wahr et al. (1981) UPDATE HISTORY: + Updated 11/2024: moved love number calculator to arguments Written 09/2024 """ import pytest import numpy as np import pyTMD.arguments -import pyTMD.predict def test_love_numbers(): """ @@ -46,7 +46,7 @@ def test_love_numbers(): for c, v in exp.items(): # calculate Love numbers omega, = pyTMD.arguments.frequency(c) - h2, k2, l2 = pyTMD.predict._body_tide_love_numbers( + h2, k2, l2 = pyTMD.arguments._love_numbers( omega, model='1066A') # check Love numbers assert np.isclose(h2, v[0], atol=15e-4)