Coverage for src/swiift/api/spectra.py: 100%
87 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-09-11 16:30 +0200
1"""Consistant implementations of usual spectral formulations.
3Spectra are defined as function of the frequency. Definitions are taken
4from [1]_.
6References
7----------
8.. [1] Stansberg, Carl Trygve & Contento, G. & Hong, S.W. & Irani, M. &
9 Ishida, S. & Mercier, R.. (2002). The specialist committee on waves
10 final report and recommendations to the 23rd ITTC. Proceedings of the
11 23rd ITTC. 505-736.
13"""
15import abc
16import functools
17import typing
19import attrs
20import numpy as np
21import scipy.integrate as integrate
23from ..lib.constants import PI_2
24from ..lib.physics import _demote_to_scalar
26"""Constant appearing in the definition of the PM family spectra."""
27_ALPHA = 8.1e-3
28"""JONSWAP peak width, below and beyond peak frequency."""
29_PEAK_WIDTH_LEQ = 0.07
30_PEAK_WIDTH_GT = 0.09
32T = typing.TypeVar("T", float, np.ndarray[tuple[int], np.dtype[np.floating]])
35class _ParametrisedSpectrum(abc.ABC):
36 @abc.abstractmethod
37 def density(self, frequencies: T) -> T:
38 """Spectral density as a function of frequency.
40 Parameters
41 ----------
42 frequencies : array_like
43 Frequencies at which to evaluate the density.
45 Returns
46 -------
47 density: float or np.ndarray
48 The density in m^2 s.
50 """
52 def discrete_energy(self, frequencies: np.ndarray) -> float:
53 r"""Spectral density integrated over the input frequencies.
55 The full spectral energy is obtained by integrating the
56 density over the positive half-line: this method computes
57 energy for the frequency band spanned by the input array.
59 The spectral density is computed for the input frequencies,
60 and integrated over these same frequencies using the
61 trapezoid rule, so that for an input of size :math:`N+1`:
63 .. math::
65 E(f_0, f_{N-1}) \approx \frac{1}{2} \sum_j=0^{N-1} {
66 (S(f_j) + S(f_{j+1})) (f_{j+1} - f_j)
67 }.
70 It is assumed the input array is sorted.
72 Parameters
73 ----------
74 frequencies : np.ndarray
75 Frequencies to use for computing the energy.
77 Returns
78 -------
79 float
80 Wave spectral energy in m^2.
82 """
83 return np.trapezoid(self.density(frequencies), frequencies) 1b
85 def _density_ang(self, angular_frequency: T) -> T:
86 # Helper method to use angular frequency as input to the density.
87 return self.density(angular_frequency / PI_2) / PI_2 1b
89 def __call__(self, frequencies: T) -> T:
90 """Wrapper for `density`.
92 Parameters
93 ----------
94 frequencies : T
95 Frequencies in Hz.
97 Returns
98 -------
99 T
100 [TODO:description]
102 """
103 return self.density(frequencies) 1bd
106@attrs.define
107class _UnimodalSpectrum(_ParametrisedSpectrum):
108 @functools.cached_property
109 def peak_period(self) -> float:
110 """Spectral peak (modal) period.
112 Returns
113 -------
114 float
115 Peak period, in s.
117 """
118 return 1 / self.peak_frequency 1e
120 @functools.cached_property
121 def peak_ang_frequency(self) -> float:
122 """Spectral peak angular frequency.
124 Returns
125 -------
126 float
127 Peak angular frequency, in rad s^-1.
129 """
130 return PI_2 * self.peak_frequency 1e
132 @property
133 @abc.abstractmethod
134 def peak_frequency(self) -> float:
135 """Peak frequency.
137 Returns
138 -------
139 float
140 Peak frequency, in Hz.
142 """
145@attrs.define
146class _PMFamily(_UnimodalSpectrum):
147 r"""Base class for spectra of the Bretschneider family.
149 These spectra are of the form
151 .. math:: S(f) = \frac{A}{f^5}\exp{\Bigl(-\frac{B}{f^4}\Bigr)}
153 with the definitions of the parameters :math:`A` and :math:`B`
154 identifying particular forms. For clarity, we call them
155 respectively `scale` and `exp_scale`.
157 In practice, spectra are parameterised by physical variables
158 which can be linked back to these parameters. In particular,
159 the peak frequency is given by
161 .. math:: f_p = {\Bigl(\frac{4B}{5}\Bigr)}^{1/4}
163 and the significant wave height by
165 .. math:: H_s = 2\sqrt{\frac{A}{B}}.
167 Attributes
168 ----------
169 scale : float
170 exp_scale : float
172 Methods
173 -------
174 __call__
176 """
178 scale: float
179 exp_scale: float
181 @functools.cached_property
182 def swh(self) -> float:
183 r"""Significant wave height.
185 We use the spectral definition,
187 .. math:: H_s = 4\sqrt{\int_0^{+\infty} S(f) df}.
189 Returns
190 -------
191 float
192 Significant wave height, in m.
194 """
195 return 2 * (self.scale / self.exp_scale) ** 0.5 1cbd
197 @functools.cached_property
198 def peak_frequency(self) -> float:
199 return (4 * self.exp_scale / 5) ** 0.25 1acbe
201 def density(self, frequencies):
202 return self.scale / frequencies**5 * np.exp(-self.exp_scale / frequencies**4) 1abd
205@attrs.define
206class PiersonMoskowitz(_PMFamily):
207 """Class encapsulating the Pierson--Moskowitz spectrum."""
209 @staticmethod
210 def _make_scale(gravity):
211 return _ALPHA * (gravity / PI_2**2) ** 2 1ac
213 @classmethod
214 def from_swh(cls, swh: float, gravity: float) -> typing.Self:
215 """Build a spectrum parameterised by the SWH.
217 Parameters
218 ----------
219 swh : float
220 Significant wave height, in m.
221 gravity : float
222 Acceleration of gravity, in m s**-2.
224 Returns
225 -------
226 PiersonMoskowitz
227 Initialised spectrum object.
229 """
230 scale = cls._make_scale(gravity) 1c
231 exp_scale = 4 * scale / swh**2 1c
232 return cls(scale, exp_scale) 1c
234 @classmethod
235 def from_peak_frequency(cls, peak_frequency: float, gravity: float) -> typing.Self:
236 """Build a spectrum parameterised by the peak frequency.
238 Parameters
239 ----------
240 peak_frequency : float
241 Peak frequency, in Hz.
242 gravity : float
243 Acceleration of gravity, in m s^-2.
245 Returns
246 -------
247 PiersonMoskowitz
248 Initialised spectrum object.
250 """
251 scale = cls._make_scale(gravity) 1ac
252 exp_scale = 5 / 4 * peak_frequency**4 1ac
253 return cls(scale, exp_scale) 1ac
256class Bretschneider(_PMFamily):
257 """Class encapsulating the Bretschneider spectrum."""
259 @classmethod
260 def from_peak_frequency_swh(cls, peak_frequency: float, swh: float) -> typing.Self:
261 """Build a spectrum parameterised by peak frequency and SWH.
263 Parameters
264 ----------
265 peak_frequency : float
266 Peak frequency, in Hz.
267 swh : float
268 Significant wave height, in m.
270 Returns
271 -------
272 Bretschneider
273 Initialised spectrum.
275 """
276 scale = 5 * swh**2 * peak_frequency**4 / 16 1ac
277 exp_scale = 5 / 4 * peak_frequency**4 1ac
278 return cls(scale, exp_scale) 1ac
281@attrs.define
282class JONSWAP(_UnimodalSpectrum):
283 """Class encapsulating the JONSWAP spectrum.
285 The JONSWAP spectrum modifies the PM spectrum by introducing a peak
286 enhancement factor.
288 Attributes
289 ----------
290 peakedness : float
291 _base_spectrum : PiersonMoskowitz
293 """
295 peakedness: float
296 _base_spectrum: PiersonMoskowitz
298 @property
299 def peak_frequency(self):
300 """Peak frequency.
302 Returns
303 -------
304 float
305 Peak frequency, in m.
307 """
308 return self._base_spectrum.peak_frequency 1bed
310 @functools.cached_property
311 def swh(self):
312 r"""Significant wave height.
314 The SWH is computed with a cubic approximation.
316 Returns
317 -------
318 float
319 Significant wave height, in m.
321 """
322 return ( 1d
323 self._base_spectrum.scale**0.5
324 * (
325 1.555
326 + 0.2596 * self.peakedness
327 - 0.02231 * self.peakedness**2
328 + 0.001142 * self.peakedness**3
329 )
330 / self.peak_frequency**2
331 )
333 @classmethod
334 def from_parameters(
335 cls,
336 peak_frequency: float,
337 peakedness: float,
338 gravity: float,
339 ) -> typing.Self:
340 """Build a spectrum.
342 Peak frequency and gravity are used to instantiate the
343 underlying Pierson--Moskowitz spectrum, while peakedness is
344 specific to the JONSWAP formulation.
346 Parameters
347 ----------
348 peak_frequency : float
349 Peak frequency, in Hz.
350 peakedness : float
351 Peakedness (base of the peak enhancement factor).
352 gravity : float
353 Acceleration of gravity, in m s^-2.
355 Returns
356 -------
357 JONSWAP
358 Initialised spectrum.
360 """
361 pm_spectrum = PiersonMoskowitz.from_peak_frequency(peak_frequency, gravity)
362 return cls(peakedness, pm_spectrum)
364 # Decoratar actually mandatory for integrate.quad not to error out
365 # on this method.
366 @_demote_to_scalar 1abd
367 def density(self, frequencies): 1abd
368 # The dimension of `frequencies` cannot be 0,
369 # to allow for masking/assignement
370 peak_width = np.ones_like(np.atleast_1d(frequencies)) * _PEAK_WIDTH_LEQ 1bd
371 peak_width[frequencies > self.peak_frequency] = _PEAK_WIDTH_GT 1bd
372 peak_enhancement = self.peakedness ** np.exp( 1bd
373 -((frequencies - self.peak_frequency) ** 2)
374 / (2 * peak_width**2 * self.peak_frequency**2)
375 )
376 return self._base_spectrum(frequencies) * peak_enhancement 1bd