pymarine.waves namespace

Submodules

pymarine.waves.wave_fields module

Implementation of the wave field solution of 1D and 2D waves. For the Wave1D and Wave2D class, a power spectral density distribution function \(S(k)\) can be imposed on the wave nodes such that we can simulate the propagation of waves obeying either a Jonswap or Gauss spectrum.

The theoretical background has the following contents

• Wave1DDescription of the Wave1D class

  • dft1d : Discrete Fourier Transform implementation for 1D wave spectra

  • fft1d : Fast Fourier Transform implementation for 1D waves spectra

• Wave2DDescription of the Wave2D class

  • dft2d : Discrete Fourier Transform implementation for 2D wave spectra

  • fft2d : Fast Fourier Transform implementation for 2D wave spectra

• References : List of literature used for this implementation

Wave1D: Theoretical background on one-D waves

Wave equation with DFT

Assume a spectrum \(S(k_i)\), with \(i=0, n_k\). The wave field \(\eta(\mathbf{x}, t)\) at time \(t\) along spatial positions \(x_p\) (where \(p=0, n_x\)) is constructed using a Discrete Fourier Transform (DFT) as follows

\[\eta(x_p, t) = \sum_i^{n_k} a_i \exp(j (k_i x_p - \omega_i t + \phi_i))\]

where \(j^2\equiv -1\), the amplitude \(a_i\) follows from the power spectral density \(S(\mathbf{k})\) as

\[a_i = \sqrt{2 S(k_i) dk_i},\]

\(\phi_i\) is the random phase corresponding to the wave node \(i\), and the angular frequency \(\omega_i\) is related to the wave vector via the deep water dispersion relation according to

\[\omega_i = \sqrt{g |k_i|}\]

with \(g=9.81 m/s^2\) being the gravity constant. We can rewrite the wave equation in matrix form:

\[\eta(x_p, t) = \sum_i^{n_k} M_{ip}\exp(-j \omega_i t)\]

where the matrix \(M_{ip}\) is

\[M_{ip} = \hat{A}_i \exp(j k_i x_p)\]

and the complex amplitude \(\hat{A}_i\) follows from

\[\hat{A}_{i} = a_i\exp(j\phi_i)\]

From the matrix form is can be seen that the wave equation can be calculated by a dot product multiplication of the \(n_x \times n_k\) matrix \(M_{ip}\) and the vector \(\exp(-j\omega_i t)\) of size \(n_k \times 1\). The result is a \(1 \times n_x\) vector of the wave along the x-axis \(\eta(x_p, t)\). This is how the DFT is implemented in the code of Wave1D.

Note that we have described the wave equation in terms of the power spectral density (PSD) equation in the wave vector domain \(k\), denoted here as \(S(k)\). Normally the PSD is given in the angular frequency domain \(\omega\) as \(S^\prime(\omega)\) We can obtain the version in wave vector domain using the equality

\[S(k) = S^\prime(\omega)\frac{d\omega}{dk} = S^\prime(\sqrt{g k}) \frac{g}{2 k}.\]

Using the PSD in the wave vector domain has the advantage that we only have to deal with the periodicity of the discrete Fourier transform in the spatial domain x: a wave will be period after a length \(2\pi/\Delta k\). To get rid of this periodicity is we can simply use a spatial domain which larger than this repetition length.

Wave equation with FFT

In case the number of wave numbers nodes is equal to the number of spatial nodes, and we have ensured that

\[\hat{A}(-k_i)\exp(-j\omega t) = \hat{A}^\ast(k_i)\exp(j\omega t)\]

we can use the FFT algorithm for calculating the wave equation at time \(t\). This is because the FFT uses the symmetry rule that for any real valued signal \(F(x)\), the Fourier transform \(\hat{F}(k)\) by definition has the property that \(\hat{F}(-k) = \hat{F}^\ast(k)\) (where the \(\ast\) indicates the complex conjugate). The FFT is a much more efficient way to obtain a solution for the wave equation: the calculation time for a DFT of a signal with \(N\) nodes is proportional to \(N^2\), while the calculation time of an FFT is proportional to \(N\log(N)\). The FFT is used for the Wave1D solution when the wave_construction field is set to FFT.

Wave2D: Theoretical background on two-D waves

The Wave2D Class gives a description of a 2D propagating wave with a spectral density distribution given by either a Gauss or a Jonswap spectrum. The class extents the functionality of the Wave1D: it obtains all the information related to the 1D wave in k-vector domain from the wave1D attribute field which stores a Wave1D object. For all the information additionally required for the 2D description (like spreading function), extra attributes are added to the class.

Again, three ways to solve the wave field from the spectral components have been implemented:

1. DFTpolar: A 2D wave field is constructed from the multiplication of the Spectral density S(k) and the directional distribution D(theta). A Discrete Fourier Transform (DFT) is used to calculate the wave field. This is a straightforward implementation, but again, DFT is slow to calculate, especially in 2D, so only use this for testing. You can speed up the calculation by selecting wave node with an amplitude larger than a certain threshold.

2. FFT: if symmetry in spectral k domain is imposed, we can again use the FFT to calculate the wave field. Recommended

3. DFTcartesian: The exact same symmetric spectral amplitudes as the FFT is used, but the wave is calculated with a DFT. This is so slow that it is not possible to include all wave components, so only used for testing purposes.

2D DFT implementation

A 2D wave can be described by its spectral density distribution \(S(k)\) and the directional distribution \(D(\theta)\) as

\[E_{k, \theta}(k, \theta) = S(k)D(\theta)\]

where for \(S(k)\) we can either take a Jonswap or a Gauss distribution as implemented in the wave_spectra module in spectrum_jonswap and spectrum_gauss. Since we are working in the wave vector domain \(k\), we must first transform the spectrum from \(\omega\) to \(k\) space, just as described in the 1D wave theory. For directional distribution, the spreading_function is recommended.

The DFTpolar implementation constructs the complex amplitudes by multiplying the spectral density \(E_{k, \theta}\) with its bin width \(dk d\theta\):

\[\hat{A}_{k, \theta} = \sqrt{2 E_{k, \theta}(k, \theta) dk d\theta} \exp({j\phi})\]

where the phase \(\phi\) is a random number between \(0\sim 2\pi\). The wave amplitude is then calculate by simply summing over all the wave components we have added, similar to our 1D version

\[\eta(\mathbf{x}, t) = \sum_p^{n_{k}}\sum_q^{n_{\theta}} \hat{A}_{k_p, \theta_q} \exp(j (\mathbf{k}_{p,q}\cdot\mathbf{x}- \omega(k_p) t))\]

where \(\mathbf{k} = (k_x, k_y)\) is the 2D wave vector in cartesian coordinates belonging to the sample wave vector in polar coordinates \((k_p, \theta_q)\) with \(k=|\mathbf{k}|\). The wave vector magnitude, \(\mathbf{x}\) is the cartesian vector in spatial domain. Also we are using the deep water dispersion relation again to calculate the angular wave frequency \(\omega\) for a given wave vector \(k\).

The DFT implementation can be sped up significantly by smartly selecting nodes, for instance based on the components above a certain threshold, or we could make a non-uniform distribution. Both methods have been implemented and can be selected by settings the wave_selection attribute field to Subrange or EqualEnergyBins, respectively.

2D FFT implementation

As already mentioned, the DFT implementation is very slow. We can speed up by selecting wave nodes in a smart way, but this consequently also takes out information from our resulting wave field. On top of that: it turns out that even after selecting only the most important spectral wave components, the DFT still is not very fast.

The recommended way to construct a wave field is again by using the FFT setting for the wave_construction argument. There is only one catch: the FFT algorithm requires the spectral components described in the cartesian wave vector domain, whereas so far we have been working in the polar wave vector domain. On top of that, we have to impose the symmetry required for a FFT:

\[\hat{A}(-\mathbf{k})\exp(-j\omega t)= \hat{A}^\ast(\mathbf{k})\exp(j\omega t)\]

This symmetry rule is the 2D version from the one we have seen in the 1D wave section.

For the FFT method we need a description of the complex amplitudes in the Cartesian domain, i.e. we need to derive \(E_{k}(\mathbf{k})\) from our starting point E_{k, theta}(k, theta). It can be shown that we need to do

\[E_{k}(\mathbf{k}) = \frac{E_{k, \theta}(k, \theta)}{k}\]

This equation is used in the function spectrum2d_complex_amplitudes as we divide the directional spreading function with the wave vector magnitude. For a more detailed derivation, please have a look at the References.

References

• Det Norske Veritas: Environmental conditions and environmental loads, April 2007, DNV-RP-C205

• Jocely Fr’echot: Realistic simulation of ocean surface using wave spectra. Proceedings of the First International Conference on Computer Graphics Theory and Applications (GRAPP, 2006), 2006, Portugal, p 76–83 <hal-00307938>

Start of the wave_fields implementation

class pymarine.waves.wave_fields.PlotProperties(show=True, color='w', linewidth=1, scattersize=0)

Bases: object

A class containing the properties of the line plot.

Parameters:

• show (bool, optional,) – Show the current wae

• color (str, optional) – Store the color of the line plotting the wave. Default = “w”

• linewidth (int, optional) – Width of the line plotting the wave. Default = 1

• scattersize (int, optional) – Size of the scatter points plotted at the wave. Default = 0, i.e., no scatter points are plotted

Notes

• Only used by WaveSimulator to connect plot properties to each wave via the plot attribute field of the wave

• Not critical and perhaps this class should be moved out of the wave_fields module later. For now, leave it as it is. It doesn’t hurt me as well.

class pymarine.waves.wave_fields.Wave1D(name=None, xmin=0, Lx=None, xmax=None, t_start=0, delta_t=1, t_length=None, nt_samples=1000000, nx_points=512, E_limit_low=0.01, E_limit_high=0.9, kx_min=0.0, kx_max=1.5708, delta_kx=0.06283, n_kx_nodes=512, n_bins_equal_energy=64, use_subrange_energy_limits=False, sample_every=1, Hs=3.0, Tp=10.0, gamma=3.3, sigma=0.0625, random_phase=False, wave_construction='FFT', wave_selection='All', lock_nodes_to_wave_one=False, spectrum_type='jonswap', spectral_version='sim', gravity0=9.80665)

Bases: object

A class for linearized solutions of 1D wave vector equation

\[\eta(x_p, t) = \sum_i^{n_k} a_i \exp(j (k_i x_p - \omega_i t + \phi_i))\]

where the amplitudes \(a_i\) follow from the Jonswap or Gauss power spectral density function. See the top of the module for more information

Parameters:

• xmin (float, optional) – Start of spatial domain. Default = 0

• Lx (float, optional) – Length of spatial domain. Default = 100

• t_start (float, optional) – Start of temporal domain. Default = 0

• t_length (float, optional) – Length of temporal domain. Default = 100

• nt_samples (int, optional) – Number of time samples. Default = 256

• nx_points (float, optional) – Number of space domain samples in x direction. Default = 256

• E_limit_low (float, optional) – Fraction of energy to cut at the low end of the spectrum. Default = 0.01

• E_limit_high (float, optional) – Fraction of energy to keep at the high end of the spectrum. Default = 0.9 (so 0.1 is cut off at the high end)

• kx_min (float, optional) – Minimal value of the wave vector domain. Default = 0

• kx_max (float, optional) – Maximum value of the wave vector domain. Default = pi/2

• delta_kx (float, optional) – Spacing in wave vector domain. Default = 0.06283 ((kx_max - kx_min) / n_kx_nodes)

• n_kx_nodes (int, optional) – Number of nodes in the wave vector domain. Default = 256

• n_bins_equal_energy (int, optional) – Number of nodes in the wave vector domain in case we are in the EqualEnergyBins mode. See wave_selection below. Default value = 100

• Hs (float, optional) – Significant wave height [m]. Default = 3 m

• Tp (float, optional) – Peak period [s]. Default = 10 s

• gamma (float, optional) – Peakness factor in case Jonswap spectrum is used. Default = 3.3

• sigma (float, optional) – Width of spectrum in case Gauss spectrum with spectral_version=’dnv’ is used. Default = 0.0625

• random_phase (bool, optional) – Use random phase with each run. Default = False

• spectrum_type ({"jonswap", "gauss"}, optional) – Type of spectrum used. Default = “jonswap”

• spectral_version ({"sim", "dnv"}, optional) – Version of wave spectrum used. Default = “sim”. For the Jonswap spectrum there is virtually no difference between the “sim” and “dnv” version. For the Gauss spectrum the “sim” version has a spectral width related to Tp which is fixed, wheras the “dnv” version has a width based on the sigma input argument.

• gravity0 (float, optional) – Gravitation constant. Default = g0 = 9.81

• wave_construction ({"FFT", "DFTpolar", "DFTcartesian"}) –

  Method how the wave field is constructed from the spectrum. Default = “FFT”. The options are

  • DFTpolar: A DFT (see top of module) is used to calculate the wave field. Slow, however, the location of the wave vectors is not restricted. Therefore, a smart selection of the wave nodes can be made. See wave_selection for the possibilities.

  • FFT : The wave is constructed from the wave vector nodes using a Fast Fourier Transfor (FFT). This implies that the spectrum nodal poits are restricted: the number of wave nodes equals the number of spatial points x, and more over, the complex amplitude must obey the symmetry rule \(a(k)=a^*(-k)\). The implies that we can not make a wave selection. However, the FFT algorithm is so much faster than the DFT that we can still obtain faster calculations time using the FFT based on all wave nodes. The FFT wave construction is therefore recommended.

  • DFTcartesian: This choice is for validation purpose only. It assumes the same symmetric spectrum as used for the FFT option but then still the slow DFT is used to calculate the wave field.

• wave_selection ({"All", "EqualEnergyBins", "Subrange"}) –

  For the DFTPolar wave construction mode we can make a selection of wave components in order to speed up the wave calculation. Three choices are possible

  • All: No selection is made so all the wave vectors as defined in the kx_nodes domain are used

  • Subrange: A Subrange of wave vectors based on the E_limit_low and E_limit_high values is made

  • EqualEnergyBins: The energy per bin is assumed equal

• use_subrange_energy_limits (bool) – If true the wave selection in the EqualEnergyBins mode is also limited by the Subrange settings. Default = True

• sample_every (int) – Make a wave selection by taking every ‘sample_every’ point in the wave vector domain. Only applicable when the wave_selection modes is Subrange

kx_nodes

Wave vector nodes of size \(n_k\) used to build the spectrum and the wave field. The values depend on the wave_selection mode and wave_construction type. In case FFT is used, the kx_nodes are symmetrical around 0 to be able to define the symmetric spectrum.

Type:

ndarray

xpoints

Spatial position array of size \(n_x\)

Type:

ndarray

amplitude

Wave height along x-direction at time t of size \(n_x\)

Type:

ndarray

phase

Random phase array of size \(n_k\)

Type:

ndarray

omega_dispersion

Angular frequency per wave node of size \(n_k\) following from the deep water dispersion relation

Type:

ndarray

spectrumK

Power spectral density of the wave spectrum along the wave vector nodes k of size \(n_k\)

Type:

ndarray

spectrumW

Power spectral density of the wave spectrum along the angular frequencies of size \(n_k\). Note that since the angular frequency are coupled to the wave vector nodes via the dispesion relation, the bin spacing of this vector is not equidistant

Type:

ndarray

complex_amplitudes

The complex amplitudes following from the power spectral density spectrumK and the phase vector as S_k * exp(j phase)

Type:

ndarray

delta_x

Spatial resolution of the x-domain

Type:

float

kx_nyquist

Maximum wave vector following from the spatial resolution delta_x as pi / delta_x

Type:

float

Examples

A wave can be setup a wave with all the default settings and plot the spectrum and wave at time 0 as follows

>>> wave_dft = Wave1D()
>>> fig, axis = wave_dft.plot_spectrum()
>>> fig, axis = wave_dft.plot_wave()

This will create the wave based on the DFT, hence it is slow. A wave construction based on the FFT can be picked as

>>> wave_fft = Wave1D()

An animation of the wave can be made as follows

>>> movie = wave_fft.animate_wave()

This will infinitely simulate the wave. We can put and finit time by setting the t_length to the wave

>>> wave_fft.reset_time(t_length=100)

For more examples, please have a look at the notebook linked via the Example section example_wave_fields

animate_wave(x_min=None, x_max=None, y_min=None, y_max=None, figsize=(12, 4), interval=25, x_time_label=0.01, y_time_label=0.93, plot_title=None, x_plot_title=0.5, y_plot_title=0.97, title_font_size=12, title_horizontal_algnment='center')

Animate the 1D wave vs space and wave vs time

Parameters:

• fig (figure) – Figure returned by previous call can be update

• ax (Axes) – Axes returned by previous call can be updated to add more waves

• x_min (float, optional) – Start of the x-range. Default = None

• x_max (float, optional) – End of the x-range. Default = None

• y_min (float, optional) – Start of the y-range. Default = None

• y_max (float, optional) – End of the y-range. Default = None

• figsize (tuple, optional) – Tuple with (with, height) in inches. Default = None

• interval (int, optional) – Time in ms between each frame. Default = 25

• x_time_label (float, optional) – X position relative to the current graph where to put the time label

• y_time_label (float, optional) – Y Position relative to the current graph where to put the time label

Returns:

(fig, ax) reference to the figure and its axis

Return type:

tuple

calculate_omega_dispersion()

Calculate the omega frequency belonging to the wave vectors according to the deep water dispersion relation.

calculate_spectra_modulus()

This routine calculates the Spectrum in omega and k domain. The spectrum can be either a Jonswap or a Gauss spectrum depending on the wave_type argument Based on the energy limits, the minimum and maximum frequency/wave vector is obtained using the specspecs routine. Based on those limits, you may take a selection of wave vectors if the wave_selection variable is set to ‘subrange’ or ‘EqualEnergyBins’

calculate_wave_surface()

Calculate the wave surface for current time using either DFT or FFT

Notes

The wave_construction attribute determines which algorithm is used to calculate the the wave surface. The following options are possible

• DFTpolar: a DFT (Discrete Fourier Transform) is used to calculate the wave surface. There are no requirement to the number and shape of the fourier nodes stored in the ‘complex_amplitudes’.

• FFT: a FFT (Fast Fourier Transform) is used to calculate the wave surface. This implies that the spectral nodes must be symmetrical around zero as S(k)=S^*(-k)

• DFTcartesian: A DFT is used to calculate the wave surface, however, still the same symmetric spectrum as used with the FFT is supplied. This allows to compare the DFT and DFT in calculation time and outcome with the exact same outcome (as the input nodes can be the same). The scaling with 0.5 due to the double spectrum (as it is symmetric) is taken care of here.

static dft_complex_amplitudes(S_tilde, exp_kx, omega, time)

Calculate the wave height from the complex amplitudes at given time

Parameters:

• S_tilde (complex ndarray) – vector Nx1 with the complex amplitudes

• exp_kx (ndarray array) – NxM matrix with exp(k * x) where k is Nx1 k waves vectors and x is 1XM spatial nodes

• omega (ndarray) – N ndvector

• time (float) – Current time

Returns:

Mx1 real array with the wave amplitude for given time

Return type:

ndarray array

Notes

Calculate the vector height(x0,x1,..xM) = real(sum_k A*exp(j(k*x-w*t))), where the MxN matrix exp(k*x) and the Nx1 vector A (with the complex amplitudes) have been pre-calculated and exp(j*w*t) is put here in a NxN diagonal matrix. Return the Mx1 vector with height for each position x

static fft_amplitude(S_tilde, omega, time)

Calculate the amplitude at time using the FFT

Parameters:

• S_tilde (ndarray) – Nx1 array with the complex amplitudes of the wave spectrum

• omega (ndarray) – N x 1 array with the real angular frequency rad/s per wave vector k

• time (float) – Current time

Returns:

real Nx1 array with the wave amplitude at time t

Return type:

ndarray

Notes

Since the FFT is used, the S_tilde should by symmetrical around k=0 such that S(k)=S^*(-k)

static limit_markers(ax, text, a_x, a_y, s_x=0, s_y=0, symbol='vy', arrow_color='black', fontsize=8)

make_report()

Make report of settings for this wave

next_time()

Increase the time

plot_spectrum(fig=None, axis=None, x_min=None, x_max=None, y_min=None, y_max=None, figsize=None, add_limit_markers=False, add_hs_estimate=False, add_n_points=True, y_n_points_label=0.8, line_markers='o', markersize=4, linestyle='-', linecolor=None, plot_title=None, x_plot_title=0.5, y_plot_title=0.97, title_font_size=12, title_horizontal_algnment='center')

Plot the spectrum in k-spaco

Parameters:

• fig (Figure, optional) – Reference to the figure to add the plot to

• axis (Axis, optional) – Reference to the axis

• x_min (float, optional) – set minimum x value, Default = None i.e take the default

• x_max (float, optional) – set maximum x value, Default = None i.e take the default

• y_min (float, optional) – set minimum y value, Default = None i.e take the default

• y_max (float, optional) – set maximum y value, Default = None i.e take the default

• figsize (tuple) – (width, height) of Figure. Default = None ie. take the default

• add_limit_markers (bool, optional) – Add markers to indicate the low and high limits of the energy in the spectrum

• add_hs_estimate=False – Add a label indicating the Hs belonging to the energy of the spectrum

• optional – Add a label indicating the Hs belonging to the energy of the spectrum

• add_n_points (bool, optional) – Add a label with the number of points

• y_n_points_label (bool, optional) – Give the y coordinates the n_point label, Default = 0.8

• line_markers (str, optional) – Give the marker of the plot. Default = “o”. Set to “” to remove the markers

• linestyle (str, optional) – Give the linestyle, Default = “-”

• linecolor (str, optional) – Give the color of the line. Default = None, ie take the default

• markersize (int, optional) – Give the size of the markers. Default = 4

Returns:

(fig, axis), reference to the create fig and axis

Return type:

tuple

plot_wave(fig=None, ax=None, x_min=None, x_max=None, y_min=None, y_max=None, figsize=(12, 4), linestyle='-', line_markers=None, linecolor=None, markersize=5, label=None, plot_title=None, x_plot_title=0.5, y_plot_title=0.97, title_font_size=12, title_horizontal_algnment='center')

Plot the 1D wave vs space and wave vs time

Parameters:

• fig (figure) – Figure returned by previous call can be update

• ax (Axes) – Axes returned by previous call can be updated to add more waves

• x_min (float, optional) – Start of the x-range. Default = None

• x_max (float, optional) – End of the x-range. Default = None

• y_min (float, optional) – Start of the y-range. Default = None

• y_max (float, optional) – End of the y-range. Default = None

• figsize (tuple, optional) – Tuple with (with, height) in inches. Default = None

• plot_title (str, optional) – Add a title to the plot. Default = None, ie. don’t add a title

• x_plot_title (float, optional) – X Location of the plot title. Default = 0.5 wro the figurek

• y_plot_title (float, optional) – Y Location of the plot title. Default = 0.97 wro the figurek

• title_font_size (int, optional) – Font size of the title

Returns:

(fig, ax) reference to the figure and its axis

Return type:

tuple

propagate_wave()

Increase the time and recalculate the surface

reset_time(t_length=None, t_start=0, nt_samples=10000000, delta_t=1)

Reset all time properties and allow to recalculate

set_wave_construction(mode)

Set the wave how the wave field is constructed from the power spectral density: via FFT or DFT

Parameters:

mode ({"FFT", "DFTpolar", "DFTcartesian"}) – Construction type of the wave field. In case FFT or DFTcartesian is chosen, the spectrum most be symmetric and therefore mirror must be True

update_x_k_t_sample_space()

After a change of number of x-points or wave vector nodes k has been made, call this routine to update all the arrays

class pymarine.waves.wave_fields.Wave2D(wave1D, name=None, Lx=200, Ly=200, xmin=0, ymin=0, nx_points=64, ny_points=64, kx_min=0.0, kx_max=3.1415, delta_kx=0.0981, n_kx_nodes=32, ky_min=0.0, ky_max=3.1415, delta_ky=0.0981, n_ky_nodes=32, theta_min=0, theta_max=6.283185307179586, n_theta_nodes=100, Theta_0=0, Theta_s_spreading_factor=5)

Bases: object

A class for linearized solutions of the 2D wave (deep water, linear). The Wave1D is used to describe the radial direction

Parameters:

• wave1D (Wave1D) – The Wave1D field to describe the spectral properties along the radial direction in k-domain the 2D wave

• Lx (float, optional) – Length of the domain in x-direction in m, default = 200 m

• Ly (float, optional) – Length of the domain in y-direction in m, default = 200 m

• xmin (float, optional) – Starting coordinate of the x-domain. Default = 0 m

• ymin (float, optional) – Starting coordinate of the y-domain. Default = 0 m

• nx_points (int, optional) – Number of nodes in the x-domain. Default = 64

• ny_points (int, optional) – Number of nodes in the y-domain. Default = 64

• kx_min (float, optional) – Start of wave vector domain. Default = 0 rad/m

• kx_max (float, optional) – End of wave vector domain. Default = pi rad/m

• delta_kx (float, optional) – Bin spacing in wave vector domain. Default is pi/32

• n_kx_nodes (int, optional) – Number of bins in in x direction in k vector domain. Default = 64

• ky_min (float, optional) – Start of wave vector domain. Default = 0 rad/m

• ky_max (float, optional) – End of wave vector domain. Default = pi rad/m

• delta_ky (float, optional) – Bin spacing in wave vector domain. Default is pi/32

• n_ky_nodes (int, optional) – Number of bins in y direction in k vector domain. Default = 64

• theta_min (float, optional) – Minimum value of theta direction in polar domain. Default = 0 rad

• theta_max (float, optional) – Maximum value of theta direction in polar domain. Default = 2 pi rad

• n_theta_nodes (int, optional) – Number of nodes in theta direction. Default = 100

• Theta_0 (float, optional) – Main wave direction in rad. Default = 0 rad. The valid range is: -pi/2 < T < pi

• Theta_s_spreading_factor (float, optional) – Spreading factor (s-definition) in theta domain. Default = 5 (Typical for wind waves, for swell use 13.0)

animate_wave(figsize=None, plot_title=None, x_plot_title=0.05, y_plot_title=0.99, min_data_value=None, max_data_value=None, number_of_contour_levels=10, title_font_size=10, x_time_label=0.05, y_time_label=0.95, x_hs_label=0.05, y_hs_label=0.92, add_hs_estimate=True, color_map=<matplotlib.colors.LinearSegmentedColormap object>, zorder=0, use_contourf=True, interval=1, title_horizontal_algnment='center')

Create a plot of the current wave

Parameters:

• figsize (tuple) – x and y size of the total figure. Default is None, so figure size is not imposed

• plot_title (str) – Title to put in th figure. If None, use the titles found in the liftdyn file

• x_plot_title (float, optional) – x position (as a fraction of the sub plot index_title_axis). Default = 0.0

• y_plot_title (float, optional) – y position (as a fraction of the sub plot index_title_axis). Default = 1.1

• delta_y_titles (float, optional) – spacing between the lines of thet title lines. Default = 0.05

• max_title_lines (int or None) – Maxinmum number of title lines to plot. Default = None, ie. plot all

• title_font_size (int) – Size of the title font

• r_axis_lim (tuple or None) – The limits of the x axis. Default = None, so no limits are imposed

• zorder (int, optional) – Position of the contour plot. Default = 0, meaning that the contours are placed at the bottom and the grid and labels are visible

• use_contourf (bool, optional) – If true, use contourf to make the contour plot. Slower. Default = False

Returns:

Handle to the figure and the axis

Return type:

tuple (fig, axis)

calculate_omega_dispersion()

calculate_spectral_components()

Calculate the 2D wave density function from the current k-array and spreading angles

calculate_spreading_function()

Calculate the spreading function

calculate_wave_surface()

dft_complex_amplitudes(S_tilde, omega, time)

Calculate DFT of complex amplitudes at time ‘time’

Parameters:

• S_tilde (Complex amplitude obtained from)

• omega (angular frequency for each node)

• time (current time)

Returns:

2D array with Discrete fourier transform

Return type:

ndarray

Notes

• The trick with the exponential matrix exp (j*x*k) does not work in 2D because you run out of memory too fast. Therefore, calculate the wave field with a loop over the wave vectors.

• This algorithm is really slow, so you should use FFT for 2D waves!

export_complex_amplitudes(filename, exportAsHD5=True)

Export the calculated complex amplitudes to HDF 5 file

Parameters:

• filename (str) – Name of the file

• exportAsHD5 (bool, optional) – Export as HD5. Default = True. If false, the complex amplitudes are written to Ascii

fft_amplitude(S_tilde, omega, time)

Calculate Fourier transform of S using the FFT :param S_tilde: Complex array of the fouriern components :type S_tilde: ndarray, complex :param omega: Angular freqyencies belong to the wave vectors :type omega: ndarray :param time: Current time :type time: float

Returns:

real array with the DFT of the complex amplitudes

Return type:

ndarray

make_report()

Make report of settings for this wave

plot_spectrum(figsize=(12, 6), r_axis_type='wave_number', plot_title=None, x_plot_title=0.05, y_plot_title=0.99, min_data_value=0, max_data_value=None, number_of_contour_levels=10, title_font_size=10, r_axis_lim=None, polar_projection=False, shift_origin=True, color_map=<matplotlib.colors.LinearSegmentedColormap object>, color_map_phase=<matplotlib.colors.LinearSegmentedColormap object>, zorder=0, use_contourf=False, r_label_position=180, kx_min=None, kx_max=None, ky_min=None, ky_max=None)

Create a polar plot of the current spectral data

Parameters:

• figsize (tuple) – x and y size of the total figure. Default is None, so figure size is not imposed

• r_axis_type ({"wave_number", "wave_length"}, optional) – quantity to use along the radial axis. Either wave_number [rad/m] or wave_length [m]. Default = “wave_number”

• plot_title (str) – Title to put in th figure. If None, don’t add a title

• x_plot_title (float, optional) – x position (as a fraction of the sub plot index_title_axis). Default = 0.0

• y_plot_title (float, optional) – y position (as a fraction of the sub plot index_title_axis). Default = 1.1

• title_font_size (int) – Size of the title font

• r_axis_lim (tuple or None) – The limits of the x axis. Default = None, so no limits are imposed

• zorder (int, optional) – Position of the contour plot. Default = 0, meaning that the contours are placed at the bottom and the grid and labels are visible

• r_label_position (float, optional) – Angle along which we position the labels of the r-axis

• use_contourf (bool, optional) – Use the slower contourf to create a contour plot. Default = False, i.e., imshow is used

• kx_min (float, optional) – Minimum wave vector node in x direction for cartesian plot. Default = None

• kx_max (float, optional) – Maximum wave vector node in x direction for cartesian plot. Default = None

• ky_min (float, optional) – Minimum wave vector node in y direction for cartesian plot. Default = None

• ky_max (float, optional) – Maximum wave vector node in y direction for cartesian plot. Default = None

Returns:

Handle to the figure and the axis

Return type:

tuple (fig, axis)

plot_wave(figsize=None, r_axis_type='frequency', plot_title=None, x_plot_title=0.05, y_plot_title=0.99, min_data_value=None, max_data_value=None, number_of_contour_levels=10, title_font_size=10, x_time_label=0.05, y_time_label=0.95, x_hs_label=0.05, y_hs_label=0.92, add_hs_estimate=True, color_map=<matplotlib.colors.LinearSegmentedColormap object>, zorder=0, use_contourf=True)

Create a plot of the current wave

Parameters:

• figsize (tuple) – x and y size of the total figure. Default is None, so figure size is not imposed

• r_axis_type ({"frequency", "period"}, optional) – quantity to use along the radial axis. Either frequency [Hz] or period [s]. Default = “frequency”

• plot_title (str) – Title to put in th figure. If None, use the titles found in the liftdyn file

• x_plot_title (float, optional) – x position (as a fraction of the sub plot index_title_axis). Default = 0.0

• y_plot_title (float, optional) – y position (as a fraction of the sub plot index_title_axis). Default = 1.1

• delta_y_titles (float, optional) – spacing between the lines of thet title lines. Default = 0.05

• max_title_lines (int or None) – Maxinmum number of title lines to plot. Default = None, ie. plot all

• title_font_size (int) – Size of the title font

• r_axis_lim (tuple or None) – The limits of the x axis. Default = None, so no limits are imposed

• zorder (int, optional) – Position of the contour plot. Default = 0, meaning that the contours are placed at the bottom and the grid and labels are visible

• use_contourf (bool, optional) – If true, use contourf to make the contour plot. Slower. Default = False

Returns:

Handle to the figure and the axis

Return type:

tuple (fig, axis)

propagate_wave()

Increase the time and recalculate the surface

update_k_polar_mesh()

Update the polar mesh and its bin area divided by k belonging to the current k/theta mesh

Notes

The bin area follows from the theta and k linear array as k x dtheta x dk However, to calculate the complex wave amplitude we need to divide D(omega, kk) by kk to go the cartesian mesh. Here, by not multiplying with kk to get the bin. Therefore, we calculate darea / kk = kk dtheta * dkk / kk = dtheta x dkk

update_x_k_theta_sample_space()

pymarine.waves.wave_spectra module

Marine Wave spectral distributions

Implementations of some standard wave spectra and distribution functions and other utilities.

pymarine.waves.wave_spectra.ag_taylor_expand(gamma)

Obtain the Ag parameter for the Jonswap model using a taylor expansion

Parameters:

gamma (float) – Peakness parameter of the Jonswap model

Returns:

Ag parameter

Return type:

float

pymarine.waves.wave_spectra.alpha_jonswap(wind_speed, fetch)

Calculate the alpha parameter for the Jonswap spectrum

Parameters:

• wind_speed (float) – wind speed at 10 m

• fetch (float) – fetch length in m

Returns:

alpha factor

Return type:

float

Notes

• The alpha factor is calculated as

\[\alpha_g = 0.076 \Big(\frac{{U_w}^2}{F g_0}\Big)^{0.22}\]

Where \(U_w\) is the wind velocity and \(F\) is the Fetch length

• Deprecated: it is better to define the Jonswap spectrum by explicitly specifying Hs and Tp

pymarine.waves.wave_spectra.d_omega_e_prime(omega, velocity)

Calculate the gradient of the relation between the encountered and true frequency

Parameters:

• omega (ndarray) – 1-D array with true frequency in rad/s

• velocity (float) – Velocity along the wave traveling direction.

Returns:

1-D array with the gradient of the relation between true and encountered frequency

Return type:

ndarray

Notes

For deep water waves, the relation between the true and encountered frequency is

\[\omega_e = \omega - k U = \omega - \frac{\omega^2 U}{g_0}\]

Where

Variable	Meaning
\(k\)	Wave vector given by deep water dispersion relation \(\omega^2/g_0\)
\(g_0\)	Gravity constant \(g_0=9.81\) m/s
\(U\)	Velocity parallel to wave vector direction. U > 0 is traveling with the wave

The gradient as calculated in this function is therefore

\[\frac{d\omega_e}{d\omega} = 1 - 2 \omega \frac{U}{g_0}\]

pymarine.waves.wave_spectra.initialize_phase(k_waves, seed=1)

Initialise a random phase is in the range \(0\sim 2\pi\) for all the nodes of the input vector k_waves

Parameters:

• k_waves – vector of wave components for which the phase need to be calculated

• seed (int, optional) – See the random generator with an integer such that every run you get the same random results for the same seed. Default value = 1. If seed = 0, no seed is used and this will result in a different result each random run

Returns:

Array of same length as k_waves with a random phase in the range \(0\sim 2\pi\)

Return type:

ndarray

pymarine.waves.wave_spectra.mask_out_of_range(kx, kmin, kmax)

Create a mask array with the values in between kmin and kmax True

Parameters:

• kx (ndarray) – Wave vector array

• kmin (float) – Minimum wave value

• kmax (float) – Maximum wave value

Returns:

array of length kx.size with booleans to mask the range

Return type:

ndarray

Notes

It is ensured that there is always an even amount of True values

pymarine.waves.wave_spectra.omega_critical(velocity, omega=None)

Get the critical frequency where the encountered frequency reaches its maximum value

Parameters:

• omega (ndarray, optional) – true frequency (Default value = None)

• velocity – the velocity of the monitor point in the direction of the wave vector k

• omega – array with true frequencies (Default value = None). If this array is passed, also the index of the location where the critical frequency is reached is returned

Returns:

(omega_c, i_c) Critical encountered frequency plus the location in the omega array of the critical encountered frequency. If the omega vector was not given, None is returned for i_c

Return type:

tuple

Notes

The encountered frequency based on the deep water dispersion relation is

\[\omega_e = \omega - k U = \omega - \frac{\omega^2 U}{g_0}\]

The critical frequency where the encountered frequency reaches it maximum is therefore located where the derivative of \(\omega_e\) is zero, i.e for:

\[\omega_c = \frac{g_0}{2 U}\]

pymarine.waves.wave_spectra.omega_deep_water(wave_number_vector)

Deep water wave dispersion relation

Parameters:

wave_number_vector (ndarray) – vector with the wave numbers in rad/m

Returns:

vector with omega for k. If k < 0, omega is also < 0

Return type:

ndarray

Notes

Implements the dispersion relation for deep water waves

\[\omega(k) = \sqrt{g_0 k}\]

where \(g_0=9.81~m/s^2\). A sign is added such that waves traveling in negative direction have also a negative angular frequency.

pymarine.waves.wave_spectra.omega_e_vs_omega(omega, u_monitor)

Calculate the relation between the encountered and true frequency

Parameters:

• omega (ndarray of float) – True frequency in rad/s

• u_monitor (float) – The velocity of the monitor point in the direction of the wave vector k in m/s

Returns:

Encountered frequency for omega

Return type:

float or ndarray

Notes

The following equation is used

\[\omega_e = \omega - k U = \omega - (\omega^2/g) U\]

where

\[k = \omega^2/g_0\]

and \(g_0=9.81~m/s^2\) is the gravity constant

pymarine.waves.wave_spectra.omega_peak_jonswap(wind_speed, fetch)

Calculate the angular peak frequency \(\omega_p\) for a given wind speed and fetch length

Parameters:

• wind_speed (float) – Wind speed at 10 m

• fetch (float) – Fetch length in m

Returns:

Angular peak frequency

Return type:

float

Notes

The following equation is implemented to calculate the angular peak frequency

\[\omega_p = 22 (\frac{g_0 ^ 2 }{U_w F}) ^{1/3}\]

with \(U_w\) the wind speed and \(F\) the fetch length

pymarine.waves.wave_spectra.omega_vs_omega_e(omega, velocity)

Calculate the relation between the true and encountered frequency

Parameters:

• omega (float) – Actual frequency

• velocity (float) – The velocity of the monitor point in the direction of the wave vector k

Returns:

(omega_1, omega_2): 2 real frequencies belonging to the same encountered frequency or omega_1, None if only one solution is available or (None, None) if no solutions exist

Return type:

tuple

Notes

The encountered frequency is given by

\[\omega_e = \omega - k U = \omega - (\omega^2/g) U = \omega (1 - \frac{\omega}{2\omega_c})\]

where we assume deep water waves such that

\[k = \omega^2/g_0\]

and \(\omega_c\) is the critical frequency defined as

\[\omega_c = g / (2 U)\]

Using the quadratic equation above, we can write the true angular frequency as a function of the encountered frequency as

\[\omega = \omega_c ( 1 \pm \sqrt{D})\]

with

\[D = 1 - \frac{2\omega_e}{\omega_c}\]

For \(\omega_e = \omega_c/2\) there is exactly one solution: \(\omega=\omega_c\). For \(\omega_e < \omega_c/2\), two solution are returned and for \(\omega_e > \omega_c/2\) no solutions exist.

pymarine.waves.wave_spectra.rayleigh_cdf(omega, sigma)

Calculate the Rayleigh cumulative distribution

Parameters:

• omega (class:numpy.ndarray) – 1-D array with the angular frequency in rad/s

• sigma (float) – A control parameter

Returns:

class – 1-D array with the cumulative distribution function

Return type:

numpy.ndarray

Examples

Create a array with some radial frequencies running from 0 to 2.5 rad/s

>>> frequencies = np.linspace(0, 2.5, 10)

Calculate the Rayleigh cumulative distribution with a width parameter sigma of 3.1

>>> rayleigh_cdf(omega=frequencies, sigma=3.1)
array([ 1.        ,  0.99899686,  0.99599345,  0.99100784,  0.98406987,
        0.97522096,  0.9645136 ,  0.95201092,  0.937786  ,  0.9219212 ])

Notes

The Rayleigh Cumulative Distribution is calculated as

\[C_R = \exp\Big[-2\Big(\frac{\omega}{4\sigma}\Big)^2\Big]\]

pymarine.waves.wave_spectra.rayleigh_pdf(omega, sigma)

The Raleigh probability density function.

Parameters:

• omega – array with the frequencies

• sigma – Standard deviation of signal sigma = sqrt(m_0), where m_0 is the zeroth moment

Returns:

Array of same size as omega with the Rayleigh distribution

Return type:

ndarray

Notes

The matlab uses the CDF, which is related

pymarine.waves.wave_spectra.rotate_fft_2d(data2d, angle=0, pivot_x=0, pivot_y=0)

Rotate a 2D array around a pivot

Parameters:

• data2d (ndarray) – A 2D array

• angle (float, optional) – Angle of rotation. default = 0, so no rotation

• pivot_x (int, optional) – pivot at x-axis, Default =0

• pivot_y (int, optional) – pivot at y-axis, Default =0

Returns:

Rotate 2d array of the same shape

Return type:

ndarray

pymarine.waves.wave_spectra.set_heading(data, directions, heading, heading_reverse=True, direction_axis=2)

Rotate the data with the heading of the RAO

Parameters:

• data (ndarray) – 3D array representing the ROA’s on a mesh of n_frequencies x n_direction (the directions are on axis=1). Note however, that during design phase data was passed as a 3D array in which the axis=0 direction gave the RAO, so the direction is really on axis=2. In case you have only one RAO, then the direction axis should be set to one

• directions (ndarray) – 1D array with the direction in the RAO

• heading (float) – Heading in degrees defined as where the bow is going to

• heading_reverse (bool, optional) – In LiftDyn a heading of 0 actually corresponds with stern going forward, so with this flag, we can swap the direction, .i.e. define a header like envy view as e.g. H=10 will become H=190 (Default value = True)

• direction_axis (int, optional) – Axis direction over which we roll the RAO. Default is 2

Returns:

the rolled RAO for the current heading

Return type:

type

pymarine.waves.wave_spectra.specspecs(frequency, amplitude, lowlim=0.01, higlim=0.9, mirror=False)

Calculate some specs from a spectrum: the peak frequency, total energy, etc.

Parameters:

• frequency (array_like) – the frequency axis

• amplitude (array_like) – the amplitude of the spectral values

• lowlim (float, optional) – Find the point below which the energy is a fraction lowlim of the total. Default = 0.01

• higlim (float, optional) – Find the point below which the energy is a fraction higlim of the total. Default = 0.9

• mirror (bool) – Assumed that we have a mirrored spectral distribution so use half of the domain only

Returns:

A tuple with the following values:

• f_min : frequency where spectrum first exceeds the threshold (0.01 of the peak value)

• f_peak : the location of the peak

• f01 : the frequency below with there is only 1% of the energy

• f96 : the frequency below with there is 96% of the energy. The interval f01-f96 contains 95 percent !

• variance : total energy below spectrum

• a_peak : the value at the peak location

Return type:

tuple

pymarine.waves.wave_spectra.spectrum2d_complex_amplitudes(kx_nodes, ky_nodes, Hs=1.0, Tp=10.0, gamma=3.3, sigma=0.0625, Theta_0=0.0, Theta_s_spread_kx=5, spectrum_type='jonswap', spectral_version='sim', mirror=True, seed=None)

Calculate the spectral complex amplitude on a cartesian wave vector mesh with a given spectral and directional distribution.

Parameters:

• kx_nodes (ndarray) – list of kx nodes. Should start at kx=0

• ky_nodes (ndarray) – list of ky nodes. Should start at ky=0

• Hs (float, optional) – significant wave height (Default value = 3.0)

• Tp (float, optional) – peak period (Default value = 10)

• gamma (float, optional) – peak enhancement factor (Default value = 3.3)

• sigma (float, optional) – Width of the spectrum. Only used when the spectrum_type=gauss spectrum is used

• theta_0 (float, optional) – mean direction theta (Default value = 0.0)

• theta_s_spreading (float, optional) – n factor to control spreading of waves (Default value = 5)

• spectrum_type ({"jonswap", "gauss"}) – type of spectrum used. Either “jonswap” or “gauss”. Default = “jonswap”

• spectral_version ({"sim", "dnv"}) – Which spectral distribution version to use. Default is “sim”. For the Jonswap spectrum, there is virtually no difference between the “sim” and “dnv” version. For the Gauss spectrum, the “sim” version has a fixed spectral width related to the peak period, whereas the “dnv” version allows to set the width of the Gauss via the sigma parameter. See the function spectrum_jonswap and spectrum_gauss from the wave_spectra module for more details.

• seed – seed for the random phase (Default value = None)

Returns:

• c_ampl: 2D ndarray with complex amplitudes

• s_ampl: 2D ndarray with sign to apply on exp(j omega t) component which should be conjugated along with the c_ampl components. Check for instance the Wave2D.calculate_omega_dispersion how to use this.

Return type:

(c_ampl, s_ampl)

Notes

• This function create a the complex spectral amplitudes on a wave vector mesh where the spectrum magnitude is based on a Jonswap of Gauss power spectral density and the phase is generated with a random generator.

• The mesh is cartesian, therefore only the x and y axis of the wave vector space have to be passed.

• The returned complex amplitudes have to be passed to the fft2d in order to calculate the wave height in the spatial domain. For that reason it is ensured that the complex amplitude are point symmetric around k=0 such that A(-k) = A^*(k) which is required to use the FFT

• The k-mesh must be defined such that the positive frequency are stored in the first half of the array and the negative frequencies are stored swapped on the second half of the array See for the details: https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftfreq.html

Examples

This example shows how to create the complex spetral amplitudes and how to use that to calculate a wave surface. First, let’s create the complex amplitude on based on a jonswap spectrum

>>> L = 1000
>>> nx = 64
>>> ny = 64
>>> dx = L / nx
>>> dy = L / ny
>>> kx_nodes = 2 * np.pi * np.fft.fftfreq(nx, dx)
>>> ky_nodes = 2 * np.pi * np.fft.fftfreq(ny, dy)
>>> ak, os = spectrum2d_complex_amplitudes(kx_nodes=kx_nodes, ky_nodes=ky_nodes,
...                                        Theta_0=0.0, Hs=2.3, Tp=12.2, gamma=3.3,
...                                        Theta_s_spread_kx=5.0, seed=1)

Use Parseval’s theorem to get the variance of the calculated Fourier components and turn it into a wave height using hs = 4 sqrt(variance). Note that we divide by 2, because the k-domain stores the amplitude mirrored over k=0 so each component is stored twice

>>> hs_estimate = 4 * np.sqrt(np.square(abs(ak)).sum()) / 2
>>> print("Hs = {:.1f} m".format(hs_estimate))
Hs = 2.3 m

Indeed we can see that the hs obtained from the complex amplitudes correspond with our target Hs

The complex amplitudes can be used to calculate the wave surface at a certain time t. To do that first you need to convert the wave vector numbers to a angular frequency using the deep water dispersion relation

>>> k_mesh = np.meshgrid(kx_nodes, ky_nodes)
>>> kk = np.sqrt(k_mesh[0] ** 2 + k_mesh[1] ** 2)
>>> g0 = 9.81
>>> omega_dispersion = np.sqrt(g0 * abs(kk)) * os

Now we have omega belonging to each wave vector we can calculate the amplitude of the surface

>>> n_nodes = int(ak.size / 2)
>>> for time in np.linspace(start=0, stop=10, num=6):
...    amplitudes = n_nodes * np.fft.ifft2(ak) * np.exp(1j * (-time * os))
...    hs_et2 = 4 * np.std(amplitudes)
...    print("time = {} s  Hs = {:.1f} m".format(time, hs_estimate))
time = 0.0 s  Hs = 2.3 m
time = 2.0 s  Hs = 2.3 m
time = 4.0 s  Hs = 2.3 m
time = 6.0 s  Hs = 2.3 m
time = 8.0 s  Hs = 2.3 m
time = 10.0 s  Hs = 2.3 m

The wave surface over the domain of 1000 x 1000 m2 is stored in amplitudes. The wave height Hs estimated with 4*variance indeed is 2.3 for all time steps

pymarine.waves.wave_spectra.spectrum2d_to_spectrum2d_encountered(spectrum_2d, frequencies, directions, velocity, debug_plot=False)

Apply a velocity shift on a 2d spectrum density. It

Parameters:

• spectrum_2d (ndarray) – 2d array with dimensions n_omega x n_directions

• frequencies (ndarray) – 2d array with frequencies in rad/s (created with meshgrid)

• directions (ndarray) – 2d array with directions in rad (created with mesh grid

• velocity (float) – velocity in the direction of the wave.

• debug_plot (bool, optional) – Pop up some debugging plots. (Default value = False)

Returns:

New spectrum with shifted distribution

Return type:

ndarray

Notes

It is assumed that the spectrum is defined with respect to the current heading of the ship (and not to the north). The velocity is therefore always into the theta=0 direction by definition. A peak in the spectrum at the theta=0 means that we have following waves (because k is into the positive x-axis and the u-velocity is in the positive x-axis direction)

pymarine.waves.wave_spectra.spectrum_complex_amplitudes_on_k_mesh(kx_nodes, ky_nodes, Hs=3, Tp=10, theta_0=0, theta_s_spreading=5.0, gamma=3.3, sigma=0.0625, spectral_version='sim', spectrum_type='jonswap', seed=None)

Calculate the complex amplitudes based on a spectrum with a random phase on a wave vector mesh

Parameters:

• kx_nodes (ndarray) – list of kx nodes. Should start at kx=0

• ky_nodes (ndarray) – list of ky nodes. Should start at ky=0

• Hs (float, optional) – significant wave height (Default value = 3.0)

• Tp (float, optional) – peak period (Default value = 10)

• gamma (float, optional) – peak enhancement factor (Default value = 3.3)

• sigma (float, optional) – Width of the spectrum. Only used when the spectrum_type=gauss spectrum is used

• theta_0 (float, optional) – mean direction theta (Default value = 0.0)

• theta_s_spreading (float, optional) – n factor to control spreading of waves (Default value = 5)

• spectrum_type ({"jonswap", "gauss"}) – type of spectrum used. Either “jonswap” or “gauss”. Default = “jonswap”

• spectral_version ({"sim", "dnv"}) – Which spectral distribution version to use. Default is “sim”. For the Jonswap spectrum, there is virtually no difference between the “sim” and “dnv” version. For the Gauss spectrum, the “sim” version has a fixed spectral width related to the peak period, whereas the “dnv” version allows to set the width of the Gauss via the sigma parameter. See the function spectrum_jonswap and spectrum_gauss from the wave_spectra module for more details.

• seed – seed for the random phase (Default value = None)

Returns:

Array with complex amplitudes on a 2D wave vector mesh

Return type:

ndarray

pymarine.waves.wave_spectra.spectrum_gauss(omega, Hs=1.0, Tp=10.0, sigma=0.0625, spectral_version='dnv')

Calculate the gauss spectral density function

Parameters:

• omega (ndarray) – Vector with angular frequencies [rad/s]

• Hs (float, optional) – Significant wave height in meter (Default value = 1.0 m)

• Tp (float, optional) – Peak period swell in seconds (Default value = 10.0 s)

• sigma (float, optional) – Width of the spectrum used for the spectral_version=”dnv” only. In the “sim” spectral version, the width is related to Tp. For Hs=1 m and Tp=10 s, the default value of sigma=0.0625 results in the same width for the “sim” and “dnv” spectral_version.

• spectral_version ({"sim", "dnv"}) –

  Type of Gauss spectrum used, Default=”sim”

  • ”sim”: Hs and Tp define energy of spectrum and peak location, respectively. The width of the gauss follows automatically from Tp. Larger Tp gives more narrow width of the spectrum

  • ”dnv”: Hs and Tp define energy of spectrum and peak location, respectively. The width of the gauss is explicitly defined by the sigma input argument and does not depend on Tp

Notes

• The benefit of the ‘sim’ spectral version is that the width does not need to be defined: it follows automatically from the Tp values (larger Tp leads to a more narrow width of the spectrum)

• The ‘dnv’ definition allows to set the width explicitly via the sigma input argument. This allows better control over the width of the spectrum yourself. It can also be used to define a spectrum which has all its energy focused into one signal frequency bin by setting sigma to a very small value such as 1e-6.

References

• The Orcina spectrum is defined here https://www.orcina.com/SoftwareProducts/OrcaFlex/Documentation/Help/Content/html/Waves,WaveSpectra.htm

Returns:

Vector with spectral densities equal in length to omega

Return type:

ndarray

pymarine.waves.wave_spectra.spectrum_jonswap(omega, Hs=1.0, Tp=10.0, gamma=3.3, spectral_version='dnv')

Calculate the Jonswap Wave Spectral Density function vs. the angular frequency

Parameters:

• omega (numpy.ndarray) – 1-D array with the radial frequencies in rad/s

• Hs (float, optional) – Significant wave height in meter, Default=1.0 m

• Tp (float, optional) – Peak period in seconds, Default=10.0 s

• gamma (float, optional) – Peaking factor, Default=3.3

• spectral_version ({"sim", "dnv"}, optional) – Type of Jonswap spectrum used, either “sim” or “dnv” can b used. Default=”sim”. Note that the difference is minor and can be neglected.

Returns:

1-D array with wave spectral density distribution

Return type:

ndarray

Notes

Definition Jonswap Spectrum

The coefficients defined in this function are

\[ \begin{align}\begin{aligned}\alpha = 5.0 / 16.0\\\beta = 5.0 / 4.0\\\omega_p = 2 \pi / T_p\\\omega < \omega_p : \sigma = 0.07\\\omega \ge \omega_p : \sigma = 0.09\end{aligned}\end{align} \]

The Pierson-Moskovitz spectrum is calculated as

\[S_P = \alpha {H_s}^2 {\omega_p}^4 \omega^{-5} \exp\Big(-\beta (\omega_p / \omega)^4\Big)\]

Extra parameters for Jonswap are

\[r_g = \exp\Big(-\Big[\frac{(\omega - \omega_p)^2}{(2 (\sigma \omega_p)^2}\Big])\Big)\]

How the scaling factor \(A_g\) is calculated depends on the “spectral_version” option

• spectral_version: “dnv”

  \[A_g = 1 - 0.287 \log(\gamma)\]

• spectral_version: “sim”

  \[A_g = \sum_{p=0}^{5} a_p \Big(\frac{\gamma - 4}{3}\Big)^p\]

The difference in value of \(A_g\) between both versions is numerically negligible

And finally the Jonswap spectrum is calculated as

\[S_J = A_g \gamma^{r_g} S_P\]

Coefficients Specification

The coefficient \(a_p\) for the \(A_g\) scaling factor in spectral_version=”sim” are given as

Index p	0	1	2	3	4	5
\(a_p\)	0.600	-0.214	0.084	-0.036	0.0385	-0.0247

References

• DNV-RP-C205 pg 33

pymarine.waves.wave_spectra.spectrum_jonswap_k_domain_2(k_waves, Hs=1.0, Tp=10.0, gamma=3.3, spectral_version='sim')

Explicit version of Jonswap spectrum in k domain obtained by substituting the derivative \(d\omega/dk\) into the Jonswap equation.

Parameters:

• k_waves – input array with the wave numbers

• Hs – significant wave height (Default value = 1.0)

• Tp – Peak period (Default value = 10.0)

• gamma – peaking factor (Default value = 3.3)

• spectral_version ({"sim", "dnv"}, optional) – Type of Jonswap spectrum used, either “sim” or “dnv” can b used. Default=”sim”. Note that the difference is minor and can be neglected.

Notes

• This spectrum is the same as the spectrum_wave_k_domain for spectrum_type=jonswap

• In this version the derivative between k/omega is explicitly calculated from the Jonswap definition

Returns:

Spectrum : array with spectral components relating the spectrum to the wave number

Return type:

ndarray

pymarine.waves.wave_spectra.spectrum_to_complex_amplitudes(omega, spectral_modulus, phase, mirror=False)

Turn a wave spectral density function into a set of complex amplitudes

Parameters:

• omega (ndarray) – vector of frequency components. Can be used for both k and f or omega

• spectral_modulus – vector with spectral modulus components belonging to the current omega

• phase – array with the phase components

• mirror – mirror the spectral components (Default value = False) such that we obey S(-k) = S*(k)

Returns:

array with the complex amplitudes

Return type:

ndarray

Notes

• If the mirror flag is True, the return complex amplitude are mirrored over the symmetry axis. This allows to feed the resulting complex amplitude to the fft routine

pymarine.waves.wave_spectra.spectrum_to_spectrum_encountered(spectrum, frequencies, velocity, debug_plot=False)

Apply a velocity shift on a 1d spectrum density to obtain the encountered

spectrum

Parameters:

• spectrum (ndarray) – 2d array with dimensions n_omega x n_directions

• frequencies (ndarray) – 2D array with frequencies in rad/s (created with meshgrid)

• velocity (float) – velocity in direction of wave.

• debug_plot (bool, optional) – create some plots for debugging (Default value = False)

Returns:

spectrum_e_int is the new spectrum with shifted distribution belonging to a uniformly distributed frequency array frequencies and info_matrix is a 3 x N matrix with three vectors: omega_e (encountered frequencies which are not uniformly distributed), omega_e_swap (encountered frequencies swapped around the critical frequency) and spectrum_e (encountered spectrum)

Return type:

tuple (spectrum_e_int, info_matrix )

Notes

• The encountered frequencies are

  \[\omega_e = \omega - U \omega^2 / g_0 = \omega - \frac{\omega^2}{2\omega_c}\]

  where \(\omega_c\) is the critical frequency. It can be seen that the distribution of encountered frequencies is not uniform. To deal with this, the returned spectrum_e is put on a uniform sampled frequency array which is the same as frequencies. This done by first converting the encountered spectrum to a cumulative energy, then resample this function and turn the cumulative energy into a a density function again.

• In order to remove the second branch of the encountered frequency, this solution is swapped over the critical frequency. This means that for \(\omega>\omega_c\), the solution is not valid but is still returned.

pymarine.waves.wave_spectra.spectrum_wave_k_domain(k_waves, Hs=1.0, Tp=10.0, gamma=3.3, sigma=0.0625, spectrum_type='jonswap', spectral_version='sim')

Calculate a wave spectrum in wave-vector domain based on its definition in frequency domain

Parameters:

• k_waves (ndarray) – Input array with the wave vector nodes in rad/m

• Hs (float, optional) – significant wave height (Default value = 1.0)

• Tp (float, optional) – Peak period (Default value = 10.0)

• gamma (float, optional) – peaking factor (Default value = 3.3)

• sigma (float, optional) – The width of the spectrum; only used for a Gauss spectrum J(spectrum_type=”gauss”) using the spectral_version=”dnv”. Default = 0.064 rad/s

• spectrum_type ({"jonswap", "gauss"}) – type of spectrum used. Either “jonswap” or “gauss”. Default = “jonswap”

• spectral_version ({"sim", "dnv"}) – type of spectrum_type used. Default = “sim” * “sim”: spectrum definition as used * “dnv”: spectrum definition by DNV.

Returns:

spectrum : array with spectral components relating the spectrum to the wave number

Return type:

type

Notes

• The spectrum in omega domain, \(S(\omega)\) can be converted to a spectrum in k domain, \(S(k)\) using

  \[S(k) = S(\omega) \frac{d\omega(k)}{dk} = S(\omega) \frac{1}{2}\sqrt{\frac{g_0}{k}}\]

  where the second step used the deep water dispersion relation.

• In this function the either a Jonswap of a Gauss spectrum can be used to convert to the k-domain

pymarine.waves.wave_spectra.spreading_function(theta, theta0=0.0, s_spreading_factor=5, n_spreading_factor=None)

Calculate the spreading function \(D_k(k,\theta)\) for polar coordinates according to the DNV definition.

Parameters:

• theta (ndarray) – The theta angle in radians running from 0 to 2pi

• theta0 (float, optional) – The main direction of the spreading function [rad] (Default value = 0.0 rad)

• s_spreading_factor (int) – The spreading factor \(s\). Can be any integer > 0. (Default value = 5).

• n_spreading_factor (int) – The spreading factor \(n\). Can be any integer > 0. (Default value = None) In case this argument is not None the s_spreading factor will be overruled as both are related as s = 2 n + 1

Returns:

The spreading function as a function of theta

Return type:

ndarray

Notes

• The spreading function as defined in DNV 3.5.8.4 is calculated:

  \[D(\theta)=\frac{\Gamma(1+n/2)}{\sqrt{\pi}\Gamma(1/2+n/2)} \cos^n(\theta-\theta_0)\]

  with \(n\) the spreading factor and \(\Gamma\) is the standard Gamma function.

• An alternative spreading factor \(s\) can be found in the DNV spreading function definition 3.5.8.6. :

  \[D(\theta)=\frac{\Gamma(s+1)}{2\sqrt{\pi}\Gamma(s+1/2)} \cos(\frac{1}{2}(\theta-\theta_0))^{(2s)}\]

  Where the spreading factor \(s\) is related to the factor \(n\) as

  \[s = 2 n + 1\]

  By default, this function excepts the s spreading factor via the argument s_spreading_factor, however, the spreading factor \(n\) may be passed instead via the n_spreading_factor argument as well. In that case, the value will overwrite the n_spreading_factor argument as defined by \(n= (s-1) / 2\)

• The recommended values for respectively \(n\) and \(s\) spreading factors is given as

Spreading factor	Wind waves	Swell
\(s\)	5 ~ 9	>= 13
\(n\)	2 ~ 4	>= 6

• To prevent an overflow of the gamma function, \(\Gamma\), the logarithm of the spreading function \(\log(D(\theta))\) is calculated and then converted to its normal value by returning the exponential of \(\log(D(\theta))\)

pymarine.waves.wave_spectra.spreading_function2(theta, theta0=0.0, s_spreading_factor=5, n_spreading_factor=None)

Calculate the spreading function \(D_k(k,\theta)\) for polar coordinates

Parameters:

• theta – frequency in rad/s

• theta0 – peak frequency in rad/s (Default value = 0.0)

• s_spreading_factor (int) – The spreading factor \(s\). Can be any integer > 0. (Default value = 5).

• n_spreading_factor (int) – The spreading factor \(n\). Can be any integer > 0. (Default value = None) In case this argument is given (not None) the s_spreading factor will be overruled, as both are related according to \(s = 2 n + 1\)

Returns:

Distribution function over the theta range

Return type:

ndarray

Notes

in which

• This spreading function definition implements the version 3.5.8.6 of DNV,

  \[D(\theta) = \frac{\Gamma(1+s)}{2\sqrt{\pi}*\Gamma(0.5+s)} \cos(\frac{\theta-\theta_0}{2})^{2s}\]

  for which the input argument spreading_factor \(s\) is given by

  \[s = 2 n + 1\]

  where \(n\) is the spreading factor as used in the spreading_function implementation.

• The recommended values for respectively \(n\) and \(s\) spreading factors is given as:

Spreading factor	Wind waves	Swell
\(s\)	5 ~ 9	>= 13
\(n\)	2 ~ 4	>= 6

• Note that the spreading_factor and spreading_factor2 are two alternative formulations given by DNV and should yield almost identical results.

• Formally the spreading_factor is only based on the n-spreading definition, and the spreading_factor2 function only based on the s-spreading factor definition. In the API of the function, however, you can pass for both functions either the n or s spreading factor, internally the spreading factor is converted to the n or s definition, respectively

pymarine.waves.wave_spectra.thetaspreadspecs(theta, Dspread, areafraction=0.99)

Calculate some specs from a spectrum: the peak frequency, total energy, etc.

Parameters:

• theta (ndarray) – Array with directions

• Dspread (ndarray) – Array with distribution function with the property that the integral is one

• areafraction – (Default value = 0.99)

Returns:

A tuple contain the following values:

• i_low : Index corresponding to theta_low

• i_high : Index corresponding to theta_high

• theta_low : Frequency where spectrum first exceeds the threshold (0.01 of the peak value)

• theta_high : Frequency where spectrum again is below the threshold (0.01 of the peak value)

• theta_peak : the location of the peak of the spreading function

• D_peak : the spreading value at the peak location theta_peak

• area : the area below the spreading function supplied. For check only, the area must be 1

Return type:

tuple (i_low, i_high, theta_low, theta_high, theta_peak, D_peak, area, mask)