# -*- coding: utf-8 -*-
"""
Relevant functions for **environmental conditions**.
.. dropdown:: Copyright (C) 2023-2026 Technical University of Denmark, R.E.G. Mounet
:color: primary
:icon: law
*This code is part of the NetSSE software.*
NetSSE 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, either version 3 of the License, or (at your
option) any later version.
NetSSE 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.
You should have received a copy of the GNU General Public License along
with this program. If not, see https://www.gnu.org/licenses/.
To credit the author, users are encouraged to use below reference:
.. code-block:: text
Mounet, R. E. G., & Nielsen, U. D. NetSSE: An open-source Python package
for network-based sea state estimation from ships, buoys, and other
observation platforms (version 2.2). Technical University of Denmark,
GitLab. February 2026. https://doi.org/10.11583/DTU.26379811.
*Last updated on 20-02-2026 by R.E.G. Mounet*
"""
import numpy as np
from math import gamma
[docs]
def JONSWAP_DNV(Tp, Hs, omega, gamma="standard", h=0):
"""Computes the JONSWAP spectrum corresponding to the input sea state
parameters.
The JONSWAP spectrum is formulated as a modification of
a Pierson-Moskowitz spectrum for a developing sea state in a fetch
limited situation.
Parameters
----------
Tp : float
Peak period [s].
Hs : float
Significant wave height [m].
omega : array_like of shape (Nfreq,)
Vector of angular frequencies [rad/s].
gamma : {'standard','DNV',float}, optional
Peak shape parameter [-]. The value can be user-provided as a float.
Alternatively, if ``'standard'`` is input, then ``gamma`` will take the
standard value of 3.3., while a value ``'DNV'`` as input leads to
following the procedure 3.5.5.5 described in DNV-RP-C205.
.. tip::
Use ``gamma = 1`` to output a standard Pierson-Moskowitz spectrum.
h : float, default 0
Water depth [m]. If ``h`` is specified as input argument, then the output
JONSWAP spectrum is corrected to account for finite water depth, becoming
a standard TMA spectrum as per Bouws et al. (1985).
Returns
-------
S_J : array_like of shape (Nfreq,)
Standard wave spectrum [m^2.s/rad].
References
----------
1. DNV-RP-C205, "Environmental Conditions and Environmental Loads,
April 2007.
2. Bouws, E., Gunther, H., Rosenthal, W., & Vincent, C. L. (1985).
*Similarity of the wind wave spectrum in finite depth water. 1. Spectral form*.
Journal of Geophysical Research-Oceans, 90(NC1), 975–986.
https://doi.org/10.1029/JC090iC01p00975
See Also
--------
lin_disprel : A fast and accurate approximation of the linear wave dispersion
relationship in finite water depth.
Example
-------
>>> S_J = JONSWAP_DNV(Tp,Hs,omega,gamma='standard',h=0)
"""
g = 9.81
omega_p = 2 * np.pi / Tp # angular spectral peak frequency [rad/s]
omega = np.array(omega)
# Pierson-Moskowitz spectrum:
S_PM = (
(5 / 16 * Hs**2 * omega_p**4)
* omega ** (-5)
* np.exp(-5 / 4 * (omega / omega_p) ** (-4))
)
# JONSWAP spectrum:
crit = Tp / np.sqrt(Hs)
if gamma == "standard":
gamma = 3.3
elif gamma == "DNV":
if crit <= 3.6:
gamma = 5
print("JONSWAP spectrum should be used with caution for the given (Tp,Hs)")
elif crit >= 5:
gamma = 1
print("JONSWAP spectrum should be used with caution for the given (Tp,Hs)")
else:
gamma = np.exp(5.75 - 1.15 * crit)
print("JONSWAP spectrum with gamma =", gamma)
A_gamma = 1 - 0.287 * np.log(gamma)
sigma_a = 0.07
sigma_b = 0.09
sigma = sigma_a * np.ones(np.shape(omega)) # spectral width parameter [n.d.]
sigma[np.where(omega > omega_p)] = sigma_b
S_J = (
A_gamma
* S_PM
* gamma ** (np.exp(-0.5 * ((omega - omega_p) / (sigma * omega_p)) ** 2))
)
if h > 0: # case of finite water depth
k0 = lin_disprel(omega, h)
phi = np.cosh(h * k0) ** 2 / (np.sinh(h * k0) ** 2 + h / g * omega**2)
S_J = S_J * phi # Finite water depth TMA spectrum
return S_J
[docs]
def lin_disprel(input, h=None, direction="om2k"):
"""Computes the wavenumbers from wave angular frequencies, or vice-versa, using the linear dispersion relationship.
If the water depth is specified, then the finite water depth is accounted for. Otherwise, deep water is assumed.
In the "om2k" direction in finite water depth, this function implements a fast and accurate approximation of the linear wave dispersion relationship. In the "k2om" direction, the function always uses the exact linear dispersion relationship.
Parameters
----------
input : array_like or float
The input, which can either be the wave angular frequencies [rad/s] or the wavenumbers [1/m], depending on the ``direction``.
h : array_like or float, optional
Mean water depth [m]
direction : {'om2k','k2om'}, optional
Direction of the dispersion relation to be solved. If ``'om2k'``, then the function
computes the wavenumbers from the angular frequencies. If ``'k2om'``, then the function
computes the angular frequencies from the wavenumbers.
Returns
-------
output : array_like or float
The corresponding wavenumbers [1/m] or wave angular frequencies [rad/s] to input, depending on the ``direction``.
Example
-------
>>> k = lin_disprel(omega, h)
"""
g = 9.81 # gravitational acceleration [m/s^2]
if h is None: # deep water case
match direction:
case "k2om":
k = np.array(input)
omega = np.sqrt(g * k)
return omega
case "om2k":
omega = np.array(input)
k = omega**2 / g
return k
case _:
raise ValueError("Invalid direction. Use 'om2k' or 'k2om'.")
else: # finite water depth case
match direction:
case "k2om":
if np.array(h).ndim == 0:
depth = h
k = np.array(input)
else:
depth = np.reshape(h, (1,) * input.ndim + h.shape)
k = np.reshape(input, input.shape + (1,) * h.ndim)
omega = np.squeeze(np.sqrt(g * k * np.tanh(k * depth)))
return omega
case "om2k":
omega = np.array(input)
if np.array(h).ndim == 0:
depth = h
omega = np.array(input)
else:
depth = np.reshape(h, (1,) * input.ndim + h.shape)
omega = np.reshape(input, input.shape + (1,) * h.ndim)
period = 2 * np.pi / omega
omega_bar = np.expand_dims(
(4 * np.pi**2 * depth) / (g * period**2), axis=omega.ndim
)
alphas = np.array([0.666, 0.445, -0.105, 0.272]).reshape(
(1,) * omega.ndim + (4,)
)
exponents = np.arange(1, 5).reshape((1,) * omega.ndim + (4,))
f = 1 + np.sum(
alphas * omega_bar**exponents,
axis=-1,
)
omega_bar = omega_bar[..., 0] # remove last dimension
lamb = np.sqrt(g * depth) * period * np.sqrt(f / (1 + omega_bar * f))
k = np.squeeze(2 * np.pi / lamb)
return k
case _:
raise ValueError("Invalid direction. Use 'om2k' or 'k2om'.")
[docs]
def wavespec1dto2d(S1d, theta, theta0, s):
"""Transforms a 1-D wave spectrum into a 2-D wave spectrum.
A cosine-2s function is used for the directional spreading distribution.
Parameters
----------
S1d : array_like of shape (Nfreq,)
1-D wave spectrum.
theta : array_like of shape (Ntheta,)
Vector of wave headings [deg] at which the 2-D spectrum must be computed.
theta0 : float
Mean wave direction [deg].
s : float
Spreading parameter [-], related to the exponent of the cosine function.
Returns
-------
S2d : array_like of shape (Nfreq,Ntheta)
2-D wave spectrum, as a function of both the frequency and the wave
heading [deg].
D_theta : array_like of shape (1,Ntheta)
Directional spreading function, as a function of the wave heading [deg].
References
----------
Naess, Arvid, and Torgeir Moan. 2010. Stochastic Dynamics of
Marine Structures. Stochastic Dynamics of Marine Structures. Vol.
9780521881555. Cambridge University Press.
https://doi.org/10.1017/CBO9781139021364.
Example
-------
>>> S2d = wavespec1dto2d(S1d,theta,theta0,s)
"""
C = (
2 ** (2 * s - 1) / np.pi * (gamma(s + 1)) ** 2 / gamma(2 * s + 1)
) # Normalizing constant
D_theta = np.reshape(C * np.cos(np.radians(theta - theta0) / 2) ** (2 * s), (1, -1))
S2d = np.reshape(S1d, (-1, 1)) * np.ones((len(S1d), len(theta)))
S2d = D_theta * S2d
return S2d, D_theta
[docs]
def fit_spec_tail(
f,
S,
fc,
exponent=5,
fmax=None,
n_tail=256,
amplitude_mode="match_fc",
transition_width=0.0,
):
"""
Extend a 1-D wave spectrum :math:`S(f)` with an :math:`f^{-n}` tail beyond :math:`f_c`.
Parameters
----------
f : array_like
1-D array of frequencies [Hz], strictly increasing.
S : array_like
1-D array of spectral densities [m^2/Hz], same length as ``f``.
fc : float
Cutoff frequency :math:`f_c` [Hz]; must be within the range of ``f``.
exponent : float
Exponent :math:`n` of the power-law tail to fit beyond :math:`f_c`. Default is ``5`` for a typical high-frequency wave spectrum tail.
fmax : float or None
Maximum frequency for the extended tail. If ``None`` (default), set to ``2fc``.
n_tail : int
Number of points to generate in :math:`(f_c, f_{\max}]`. Default is ``256``.
amplitude_mode : {"match_fc", "least_squares", "median"}
How to estimate the tail amplitude :math:`C` in :math:`S(f) \sim C f^{-n}` beyond :math:`f_c`.
Options are:
- `match_fc` (default): :math:`C = S(f_c) \cdot f_c^n` using interpolation at :math:`f_c`.
- `least_squares`: fit :math:`C` by minimizing :math:`\sum (S - C f^{-n})^2` over a band near :math:`f_c`.
- `median`: :math:`C = \mathrm{median}(S \cdot f^n)` over a band near :math:`f_c` (robust to outliers).
transition_width : float
Width of cosine blend band [Hz] below :math:`f_c` to smoothly transition
from original spectrum to the tail. If set to ``0`` (default), hard splice at :math:`f_c`.
Returns
-------
f_ext : ndarray
Extended frequency vector including original ``f`` and new tail points.
S_ext : ndarray
Spectrum values aligned with ``f_ext``.
C : float
Tail amplitude used in :math:`S_{tail} = C \cdot f^{-n}`.
Example
-------
>>> f = np.linspace(0.05, 1.0, 100)
>>> S = 0.1 * f**-3 * np.exp(- (f / 0.3)**2) # Example spectrum
>>> fc = 0.5
>>> f_ext, S_ext, C = fit_spec_tail(f, S, fc, exponent=5, amplitude_mode="least_squares", transition_width=0.1)
"""
f = np.asarray(f).astype(float)
S = np.asarray(S).astype(float)
if f.ndim != 1 or S.ndim != 1 or len(f) != len(S):
raise ValueError("f and S must be 1-D arrays of the same length.")
if not np.all(np.diff(f) > 0):
raise ValueError("f must be strictly increasing.")
if not (f.min() <= fc <= f.max()):
raise ValueError("fc must lie within the range of f.")
# Interpolate S at fc for matching or to define the band
S_fc = np.interp(fc, f, S)
# Define a band [fc - transition_width, fc] for estimating C
# If transition_width == 0, use last few points below fc as a fallback.
if transition_width > 0:
f_band_min = max(f.min(), fc - transition_width)
band_mask = (f >= f_band_min) & (f <= fc)
else:
# Use last 10% of samples below fc (at least 5 points)
below_fc = f < fc
idx = np.where(below_fc)[0]
if len(idx) < 5:
band_mask = below_fc # whatever is available
else:
k = max(5, int(0.1 * len(f))) # last 10% or at least 5 points
start = max(idx[-k], idx[0])
band_mask = np.zeros_like(f, dtype=bool)
band_mask[start : idx[-1] + 1] = True
f_band = f[band_mask]
S_band = S[band_mask]
# Compute amplitude C
if amplitude_mode == "match_fc":
C = S_fc * (fc**exponent)
else:
# Robust transform: target Y = S * f^n should be approximately constant = C
Y = S_band * (f_band**exponent)
if amplitude_mode == "least_squares":
# LS estimate of constant is simply the mean
C = float(np.mean(Y))
elif amplitude_mode == "median":
C = float(np.median(Y))
else:
raise ValueError(
"Invalid amplitude_mode. Choose 'match_fc', 'least_squares', or 'median'."
)
# Create tail frequencies
if fmax is None:
fmax = 2.0 * fc
if fmax <= fc:
raise ValueError("fmax must be greater than fc.")
f_tail = np.linspace(fc, fmax, n_tail + 1)[1:] # exclude fc to avoid duplication
# Tail spectrum
S_tail = C * (f_tail ** (-exponent))
# Smooth blending near fc if requested
if transition_width > 0:
# Apply a cosine ramp over [fc - transition_width, fc]
f_blend = f[(f >= fc - transition_width) & (f <= fc)]
S_blend_orig = np.interp(f_blend, f, S)
S_blend_tail = C * (f_blend ** (-exponent))
# Weights: 0 (orig) at fc - w → 1 (tail) at fc
w = (f_blend - (fc - transition_width)) / transition_width
w = np.clip(w, 0.0, 1.0)
w = 0.5 - 0.5 * np.cos(np.pi * w) # cosine smoothstep
S_blend = (1 - w) * S_blend_orig + w * S_blend_tail
# Replace original spectrum in the blending band
S_mod = S.copy()
mask_blend = (f >= fc - transition_width) & (f <= fc)
S_mod[mask_blend] = S_blend
# Concatenate modified original part and tail
f_ext = np.concatenate([f, f_tail])
S_ext = np.concatenate([S_mod, S_tail])
else:
# Hard splice: keep original up to fc, then tail
keep_mask = f <= fc
f_ext = np.concatenate([f[keep_mask], f_tail])
S_ext = np.concatenate([S[keep_mask], S_tail])
return f_ext, S_ext, C
