Add jonswap tail

LICENSE.md

@@ -44,3 +44,12 @@ author or the affirmer.
 
 
 [1]: https://creativecommons.org/publicdomain/zero/1.0/legalcode
+
+Libraries
+---------
+
+Libraries found in the "lib" directory are distributed under
+open source (or open source-like) licenses/agreements. Appropriate license
+agreements for each library can be found in the "lib" directory.
+
+ - [jonswap](https://github.com/haphaeu/jonswap)

stglib/core/waves.py

@@ -4,6 +4,10 @@
 import numpy as np
 import scipy.signal as spsig
 import xarray as xr
+from scipy.optimize import newton
+from tqdm import tqdm
+
+from ..lib import jonswap
 
 
 def make_waves_ds(ds):
@@ -70,6 +74,27 @@ def make_waves_ds(ds):
     spec["wp_peak"] = xr.DataArray(make_Tp(spec["pspec"]), dims="time")
     spec["k"] = xr.DataArray(k, dims=("time", "frequency"))
 
+    do_jonswap = False
+    if "jonswap" in ds.attrs:
+        if ds.attrs["jonswap"].lower() == "true":
+            do_jonswap = True
+
+    if do_jonswap:
+        print("Running statistics with JONSWAP tail")
+        thejonswaptail = [
+            make_tail_jonswap(
+                spec["frequency"], spec["Pnn"][burst, :], spec["tailind"][burst].values
+            )
+            for burst in tqdm(range(len(spec["time"])))
+        ]
+        spec["jonswap_pspec"] = xr.DataArray(thejonswaptail, dims=("time", "frequency"))
+
+        m0 = xr.DataArray(make_moment(spec["frequency"], spec["pspec"], 0), dims="time")
+        m2 = xr.DataArray(make_moment(spec["frequency"], spec["pspec"], 2), dims="time")
+        spec["jonswap_wh_4061"] = xr.DataArray(make_Hs(m0), dims="time")
+        spec["jonswap_wp_4060"] = xr.DataArray(make_Tm(m0, m2), dims="time")
+        spec["jonswap_wp_peak"] = xr.DataArray(make_Tp(spec["pspec"]), dims="time")
+
     return spec
 
 
@@ -223,6 +248,29 @@ def make_tail(f, Pnn, tailind):
         return np.hstack((Pnn[:ti], tail[ti:]))
 
 
+def make_tail_jonswap(f, Pnn, tailind):
+    """Make tail using JONSWAP spectrum"""
+    ti = tailind.astype(int)
+    tail = np.ones_like(f)
+    if np.isnan(tailind):
+        return np.ones_like(f) * np.nan
+    else:
+        Hs = 0
+        Tp = make_Tp(Pnn[:ti])
+        # Use Newton-Raphson iteration to get the best value of Hs for the fit
+        result = newton(
+            run_jonswap, 0, args=(Tp, 2 * np.pi * f, ti, Pnn), full_output=True
+        )
+        # we have the value of Hs, so run jonswap to get the spectrum
+        S = jonswap.jonswap(result[0], Tp, 2 * np.pi * f)
+        return np.hstack((Pnn[:ti], S[ti:]))
+
+
+def run_jonswap(Hs, Tp, w, ti, Pnn):
+    S = jonswap.jonswap(Hs, Tp, w)
+    return Pnn[ti] - S[ti]
+
+
 def make_mwd(freqs, dirs, dspec):
     """Create mean wave direction (EPIC 4062) variable"""
 

stglib/lib/__init__.py

@@ -0,0 +1 @@
+from .jonswap import jonswap

stglib/lib/jonswap/jonswap.py

@@ -0,0 +1,78 @@
+import numpy as np
+from numpy import pi, exp, cos
+import matplotlib.pyplot as plt
+
+
+def pierson_moskowitz(hs, tp, w):
+    """Return a Pierson-Moskowitz wave spectrum.
+    See `jonswap`.
+    """
+    # Peak frequency
+    wp = 2 * pi / tp
+
+    return 5/16 * hs**2 * wp**4 * w**-5 * exp(
+        -5/4 * (w/wp)**-4
+    )
+
+
+def jonswap(hs, tp, w):
+    """Return a JONSWAP wave spectrum.
+
+    hs: float 
+        significant wave height [m]
+    tp: float
+        spectral peak period [s]
+    w: array
+        frequency domain [rd/s]. 
+    """
+
+    # Peak frequency
+    wp = 2 * pi / tp
+
+    # sigma - spectral width parameter
+    s = np.where(w <= wp, 0.07, 0.09)
+
+    # gamma - Non-dimensional peak shape parameter
+    y = jonswap_gamma(hs, tp)
+
+    # normalising factor
+    Ay = 0.2 / (0.065 * y**0.803 + 0.135)
+
+    # Pierson-Moskowitz spectrum
+    S_pm = pierson_moskowitz(hs, tp, w)
+
+    return Ay * S_pm * y ** exp(-0.5 * ((w - wp) / s / wp)**2)
+
+
+def jonswap_is_valid(hs, tp):
+    return 3.6 < tp / hs**0.5 < 5
+
+
+def jonswap_gamma(hs, tp):
+    """Return JONSWAP gamma parameters after DNV"""
+    ratio = tp / hs**0.5
+    if ratio <= 3.6:
+        return 5.0
+    elif ratio < 5:
+        return exp(5.75 - 1.15 * ratio)
+    else:
+        return 1
+
+
+def test():
+    hs = 5.5
+    tp = 10
+    print('JONSWAP valid:', jonswap_is_valid(hs, tp))
+    Ts = np.linspace(1, 20, 1000)
+    ws = 2 * pi / Ts
+    PM = pierson_moskowitz(hs, tp, ws)
+    J = jonswap(hs, tp, ws)
+    plt.plot(Ts, J, label='jonswap')
+    plt.plot(Ts, PM, label='PM')
+    plt.legend()
+    plt.grid()
+    plt.show()
+
+
+if __name__ == '__main__':
+    test()

stglib/lib/jonswap/jonswap_license.txt

@@ -0,0 +1,12 @@
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, version 2.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    Contact the author should you have any question.
+
+    https://github.com/haphaeu/jonswap

stglib/vec/nc2waves.py

@@ -13,8 +13,9 @@ def nc_to_waves(nc_filename):
 
     spec = waves.make_waves_ds(ds)
 
-    for k in ["wp_peak", "wh_4061", "wp_4060", "pspec"]:
-        ds[k] = spec[k]
+    for k in ["wp_peak", "wh_4061", "wp_4060", "pspec", "pspec_jonswap"]:
+        if k in spec:
+            ds[k] = spec[k]
 
     dopuv = False
     if "puv" in ds.attrs: