enh: add Torsethaugen wave spectrum (#49)

Author: ali-cetin-4ss

docs/api_ref/index.rst

@@ -24,6 +24,7 @@ functions and methods.
     JONSWAP
     ModifiedPiersonMoskowitz
     OchiHubble
+    Torsethaugen
     multiply
     polar_to_complex
     RAO

docs/user_guide/concepts_utils/directional_spectrum.rst

@@ -88,4 +88,4 @@ and a quantile, ``q``.
     duration = 3 * 3600   # 3 hours
 
     mpm = spectrum.extreme(duration, q=0.37)   # most probable maximum (MPM)
-    q90 = spectrum.extreme(duration, q=0.99)   # 90-th quantile
+    q90 = spectrum.extreme(duration, q=0.90)   # 90-th quantile

docs/user_guide/standardized_spectra.rst

@@ -163,6 +163,65 @@ and :math:`q`). The total spectrum is obtained by adding together the two wave c
     spectrum.
 
 
+Torsethaugen
+------------
+The *Torsethaugen* spectrum allows you to set up a double-peaked spectrum that represents
+sea states which includes both a remotely generated swell component (with low frequency energy)
+and a locally generated wind component (with high frequency energy). The spectral spectral model
+was developed based on average measured spectra for Norwegian waters (Haltenbanken and Statfjord).
+The Torsethaugen spectrum is described by two parameter (i.e. :math:`H_s` and :math:`T_p`), and
+is given by:
+
+.. math::
+
+    S(\omega) = \alpha_{\gamma}S_1(\omega; H_{s1}, T_{p1})\gamma^{exp\left( -\frac{(\omega - \omega_{p1})^2}{2\sigma^2\omega_{p1}^2} \right)} + S_2(\omega; H_{s2}, T_{p2})
+
+where,
+
+.. math::
+
+    S_i(\omega) = \frac{A_{i}}{\omega^4} \exp(-\frac{B_{i}}{\omega^4})
+
+Furthermore,
+
+- :math:`A_i = \frac{3.26}{16} H_{si}^2 \omega_{pi}^3`
+- :math:`B_i = \omega_{pi}^4`
+- :math:`H_{si}` is the significant wave height for wave component :math:`i`.
+- :math:`\omega_{pi} = \frac{2\pi}{T_{pi}}` is the angular spectral peak frequency for wave component :math:`i`.
+- :math:`\gamma` is a peak enhancement factor.
+- :math:`\alpha_{\gamma}` is a normalizing factor.
+- :math:`\sigma` is the spectral width parameter, given by:
+
+.. math::
+    \sigma =
+    \begin{cases}
+        \sigma_a & \quad \text{if } \omega \leq \omega_{p1}\\
+        \sigma_b & \quad \text{if } \omega > \omega_{p1}
+    \end{cases}
+
+.. note::
+
+    The Torsethaugen spectrum is implemented according to the `original paper <https://www.sintef.no/globalassets/upload/fiskeri_og_havbruk/havbruksteknologi/2004-jsc-193.pdf/>`_
+    by Torsethaugen and Haver published in 2004. Refer to this paper for full implementation
+    details.
+
+The :class:`~waveresponse.Torsethaugen` class provides functionality for generating a 1-D
+Torsethaugen spectrum given two parameters (i.e., :math:`H_s`, :math:`T_p`):
+
+.. code:: python
+
+    import numpy as np
+    import waveresponse as wr
+
+
+    freq = np.arange(0.01, 1, 0.01)
+    spectrum = wr.Torsethaugen(freq, freq_hz=True)
+
+    hs = 3.5
+    tp = 10.0
+    freq, vals = spectrum(hs, tp)
+
+
 Directional spectrum
 ####################
 The directional spectrum is usually standardized in a similar way as the 1-D frequency

src/waveresponse/__init__.py

@@ -17,6 +17,7 @@
     BaseSpectrum1d,
     ModifiedPiersonMoskowitz,
     OchiHubble,
+    Torsethaugen,
 )
 from ._transform import (
     rigid_transform,
@@ -47,5 +48,6 @@
     "rigid_transform_heave",
     "rigid_transform_surge",
     "rigid_transform_sway",
+    "Torsethaugen",
     "WaveSpectrum",
 ]

src/waveresponse/_standardized1d.py

@@ -390,3 +390,203 @@ def _C(self, hs, tp, q):
     def _d(self, tp, q):
         omega_p = 2.0 * np.pi / tp
         return (4 * q + 1) * omega_p**4 / 4.0
+
+
+class Torsethaugen(BaseSpectrum1d):
+    """
+    Torsethaugen spectrum, given as:
+
+    ``S(w) = alpha * S_1(w; Hs_1, Tp_1) * gamma ** b + S_2(w; Hs_2, Tp_2)``
+
+    where,
+
+    ``S_i(w) = A_i / w ** 4 * exp(-B_i / w ** 4)``
+
+    and,
+
+    ``b = exp(-(w - wp_1) ** 2 / (2 * sigma ** 2 * wp_1 ** 2))``
+
+    Furthermore,
+
+    - ``A_i = 3.26 / 16 * Hs_i ** 2 * wp_i ** 3``
+    - ``B_i = wp_i ** 4``.
+    - ``Hs_i`` is the significant wave height of the primary and the secondary system.
+    - ``wp_i = 2pi / Tp_i`` is the angular spectral peak frequency of the primary and the secondary system.
+    - ``gamma`` is a peak enhancement factor.
+    - ``alpha = (1 + 1.1 * log(gamma) ** 1.19) / gamma`` is a normalizing factor.
+    - ``sigma`` is the spectral width parameter (established from experimental data):
+        - ``sigma = 0.07`` for ``w <= wp_1``
+        - ``sigma = 0.09`` for ``w > wp_1``
+
+    Parameters
+    ----------
+    freq : array-like
+        Sequence of frequencies to use when generating the spectrum.
+    freq_hz : bool
+        Whether the provided frequencies are in rad/s (default) or Hz.
+
+    Notes
+    -----
+    The Torsethaugen spectrum is implemented as described in reference [1]_. All the
+    intermediate parameters are applied as defined in Table 1 in reference [1]_.
+
+    References
+    ----------
+    .. [1] Torsethaugen K and Haver S, 2004. Simplified double peak spectral model
+           for ocean waves, Paper No. 2004-JSC-193, ISOPE 2004 Touson, France.
+
+    See Also
+    --------
+    ModifiedPiersonMoskowitz : Modified Pierson-Moskowitz wave spectrum.
+    JONSWAP : JONSWAP wave spectrum.
+    OchiHubble : Ochi-Hubble (three-parameter) wave spectrum.
+
+    """
+
+    # Parameters as specified in Table 1 in reference [1].
+    _af = 6.6
+    _ae = 2.0
+    _au = 25
+    _a10 = 0.7
+    _a1 = 0.5
+    _kg = 35.0
+    _b1 = 2.0
+    _a20 = 0.6
+    _a2 = 0.3
+    _a3 = 6.0
+    _g = 9.80665
+
+    def __call__(self, hs, tp, freq_hz=None):
+        """
+        Generate wave spectrum.
+
+        Parameters
+        ----------
+        hs : float
+            Significant wave height, Hs.
+        tp : float
+            Peak period, Tp.
+        freq_hz : bool, optional
+            Whether to return the frequencies and spectrum in terms of Hz (`True`)
+            or rad/s (`False`). If `None` (default), the original units of `freq` is
+            preserved.
+
+        Return
+        ------
+        freq : 1-D array
+            Frequencies corresponding to the spectrum values. Unit is set according
+            to `freq_hz`.
+        spectrum : 1-D array
+            Spectrum values. Unit is set according to `freq_hz`.
+
+        Notes
+        -----
+        The scaling between wave spectrum in terms of Hz and rad/s is defined
+        as:
+
+        ``S(f) = 2*pi*S(w)``
+
+        where ``S(f)`` and ``S(w)`` are the same spectrum but expressed
+        in terms of Hz and rad/s, respectively.
+        """
+
+        return super().__call__(hs, tp, freq_hz=freq_hz)
+
+    def _spectrum(self, omega, hs, tp):
+        tpf = self._tpf(hs)
+
+        if tp <= tpf:  # wind dominated
+            hs_primary = hs * self._rw(hs, tp)
+            tp_primary = tp
+            gamma = self._gamma(hs_primary, tp_primary)
+
+            hs_secondary = hs * np.sqrt(1.0 - self._rw(hs, tp) ** 2)
+            tp_secondary = tpf + self._b1
+
+        else:  # swell dominated
+            hs_primary = hs * self._rs(hs, tp)
+            tp_primary = tp
+            gamma = self._gamma(hs, tpf) * (1.0 + self._a3 * self._eps_u(hs, tp))
+
+            hs_secondary = hs * np.sqrt(1.0 - self._rs(hs, tp) ** 2)
+            tp_secondary = self._af * hs_secondary ** (1.0 / 3.0)
+
+        omega_p_primary = 2.0 * np.pi / tp_primary
+        A_primary = (3.26 / 16.0) * hs_primary**2 * omega_p_primary**3
+        B_primary = omega_p_primary**4
+
+        spectrum_primary = (
+            self._alpha(gamma)
+            * A_primary
+            / omega**4
+            * np.exp(-B_primary / omega**4)
+            * gamma ** self._b(tp_primary)
+        )
+
+        omega_p_secondary = 2.0 * np.pi / tp_secondary
+        A_secondary = (3.26 / 16.0) * hs_secondary**2 * omega_p_secondary**3
+        B_secondary = omega_p_secondary**4
+
+        spectrum_secondary = (
+            A_secondary / omega**4 * np.exp(-B_secondary / omega**4)
+        )
+
+        return spectrum_primary + spectrum_secondary
+
+    def _tpf(self, hs):
+        return self._af * hs ** (1.0 / 3.0)
+
+    def _tl(self, hs):
+        return self._ae * np.sqrt(hs)
+
+    def _eps_l(self, hs, tp):
+        tpf = self._tpf(hs)
+        tl = self._tl(hs)
+
+        eps_l = (tpf - tp) / (tpf - tl)
+
+        if eps_l < 0.0:
+            raise ValueError("`eps_l` not defined.")
+
+        return min(eps_l, 1.0)
+
+    def _eps_u(self, hs, tp):
+        tpf = self._tpf(hs)
+        tu = self._au
+
+        eps_u = (tp - tpf) / (tu - tpf)
+
+        if eps_u < 0.0:
+            raise ValueError("`eps_u` not defined.")
+
+        return min(eps_u, 1.0)
+
+    def _rw(self, hs, tp):
+        eps_l = self._eps_l(hs, tp)
+        return self._a10 + (1.0 - self._a10) * np.exp(-((eps_l / self._a1) ** 2.0))
+
+    def _rs(self, hs, tp):
+        eps_u = self._eps_u(hs, tp)
+        return self._a20 + (1.0 - self._a20) * np.exp(-((eps_u / self._a2) ** 2.0))
+
+    def _gamma(self, hs_, tp_):
+        s = (2.0 * np.pi / self._g) * hs_ / tp_**2.0
+        return self._kg * s ** (6.0 / 7.0)
+
+    def _alpha(self, gamma):
+        return (1.0 + 1.1 * np.log(gamma) ** 1.19) / gamma
+
+    def _b(self, tp):
+        omega_p = 2.0 * np.pi / tp
+        sigma = self._sigma(omega_p)
+        return np.exp(-0.5 * ((self._freq - omega_p) / (sigma * omega_p)) ** 2)
+
+    def _sigma(self, omega_p):
+        """
+        Spectral width.
+        """
+        arg = self._freq <= omega_p
+        sigma = np.empty_like(self._freq)
+        sigma[arg] = 0.07
+        sigma[~arg] = 0.09
+        return sigma