danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
basic
master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '7fcc4d18bd'
${ noResults }
pywafo/wafo/spectrum/models.py

2140 lines
68 KiB
Python
Raw Blame History

#!/usr/bin/env python
"""
Models module
-------------
Dispersion relation
-------------------
k2w - Translates from wave number to frequency
w2k - Translates from frequency to wave number
Model spectra
-------------
Bretschneider - Bretschneider spectral density.
Jonswap - JONSWAP spectral density
McCormick - McCormick spectral density.
OchiHubble - OchiHubble bimodal spectral density model.
Tmaspec - JONSWAP spectral density for finite water depth
Torsethaugen - Torsethaugen double peaked (swell + wind) spectrum model
Wallop - Wallop spectral density.
demospec - Loads a precreated spectrum of chosen type
jonswap_peakfact - Jonswap peakedness factor Gamma given Hm0 and Tp
jonswap_seastate - jonswap seastate from windspeed and fetch
Directional spreading functions
-------------------------------
Spreading - Directional spreading function.
"""
from __future__ import absolute_import, division
import warnings
from scipy.interpolate import interp1d
import scipy.optimize as optimize
import scipy.integrate as integrate
import scipy.special as sp
from scipy.fftpack import fft
import numpy as np
from numpy import (inf, atleast_1d, newaxis, minimum, maximum, array,
asarray, exp, log, sqrt, where, pi, arange, linspace, sin,
cos, sinh, isfinite, mod, expm1, tanh, cosh, finfo,
ones, ones_like, isnan, zeros_like, flatnonzero, sinc,
hstack, vstack, real, flipud, clip)
from ..wave_theory.dispersion_relation import w2k, k2w # @UnusedImport
from .core import SpecData1D, SpecData2D
__all__ = ['Bretschneider', 'Jonswap', 'Torsethaugen', 'Wallop', 'McCormick',
'OchiHubble', 'Tmaspec', 'jonswap_peakfact', 'jonswap_seastate',
'Spreading', 'w2k', 'k2w', 'phi1']
_EPS = finfo(float).eps
def sech(x):
return 1.0 / cosh(x)
def _gengamspec(wn, N=5, M=4):
""" Return Generalized gamma spectrum in dimensionless form
Parameters
----------
wn : arraylike
normalized frequencies, w/wp.
N : scalar
defining the decay of the high frequency part.
M : scalar
defining the spectral width around the peak.
Returns
-------
S : arraylike
spectral values, same size as wn.
The generalized gamma spectrum in non-
dimensional form is defined as:
S = G0.*wn.**(-N).*exp(-B*wn.**(-M)) for wn > 0
= 0 otherwise
where
B = N/M
C = (N-1)/M
G0 = B**C*M/gamma(C), Normalizing factor related to Bretschneider form
Note that N = 5, M = 4 corresponds to a normalized
Bretschneider spectrum.
Examples
--------
>>> import wafo.spectrum.models as wsm
>>> import numpy as np
>>> wn = np.linspace(0,4,5)
>>> wsm._gengamspec(wn, N=6, M=2)
array([ 0. , 1.16765216, 0.17309961, 0.02305179, 0.00474686])
See also
--------
Bretschneider
Jonswap,
Torsethaugen
References
----------
Torsethaugen, K. (2004)
"Simplified Double Peak Spectral Model for Ocean Waves"
In Proc. 14th ISOPE
"""
w = atleast_1d(wn)
S = zeros_like(w)
k = flatnonzero(w > 0.0)
if k.size > 0:
B = N / M
C = (N - 1.0) / M
# A = Normalizing factor related to Bretschneider form
# A = B**C*M/gamma(C)
# S[k] = A*wn[k]**(-N)*exp(-B*wn[k]**(-M))
logwn = log(w.take(k))
logA = (C * log(B) + log(M) - sp.gammaln(C))
S.put(k, exp(logA - N * logwn - B * exp(-M * logwn)))
return S
class ModelSpectrum(object):
type = 'ModelSpectrum'
def __init__(self, Hm0=7.0, Tp=11.0, **kwds):
self.Hm0 = Hm0
self.Tp = Tp
def tospecdata(self, w=None, wc=None, nw=257):
"""
Return SpecData1D object from ModelSpectrum
Parameter
---------
w : arraylike
vector of angular frequencies used in discretization of spectrum
wc : scalar
cut off frequency (default 33/Tp)
nw : int
number of frequencies
Returns
-------
S : SpecData1D object
member attributes of model spectrum are copied to S.workspace
"""
if w is None:
if wc is None:
wc = 33. / self.Tp
w = linspace(0, wc, nw)
S = SpecData1D(self.__call__(w), w)
try:
S.h = self.h
except AttributeError:
pass
S.labels.title = self.type + ' ' + S.labels.title
S.workspace = self.__dict__.copy()
return S
def chk_seastate(self):
""" Check if seastate is valid
"""
if self.Hm0 < 0:
raise ValueError('Hm0 can not be negative!')
if self.Tp <= 0:
raise ValueError('Tp must be positve!')
if self.Hm0 == 0.0:
warnings.warn('Hm0 is zero!')
self._chk_extra_param()
def _chk_extra_param(self):
pass
class Bretschneider(ModelSpectrum):
"""
Bretschneider spectral density model
Member variables
----------------
Hm0 : significant wave height (default 7 (m))
Tp : peak period (default 11 (sec))
N : scalar defining decay of high frequency part. (default 5)
M : scalar defining spectral width around the peak. (default 4)
Parameters
----------
w : array-like
angular frequencies [rad/s]
The Bretschneider spectrum is defined as
S(w) = A * G0 * wn**(-N)*exp(-N/(M*wn**M))
where
G0 = Normalizing factor related to Bretschneider form
A = (Hm0/4)**2 / wp (Normalization factor)
wn = w/wp
wp = 2*pi/Tp, angular peak frequency
This spectrum is a suitable model for fully developed sea,
i.e. a sea state where the wind has been blowing long enough over a
sufficiently open stretch of water, so that the high-frequency waves have
reached an equilibrium. In the part of the spectrum where the frequency is
greater than the peak frequency (w>wp), the energy distribution is
proportional to w**-5.
The spectrum is identical with ITTC (International Towing Tank
Conference), ISSC (International Ship and Offshore Structures Congress)
and Pierson-Moskowitz, wave spectrum given Hm0 and Tm01. It is also
identical with JONSWAP when the peakedness factor, gamma, is one.
For this spectrum, the following relations exist between the mean
period Tm01 = 2*pi*m0/m1, the peak period Tp and the mean
zero-upcrossing period Tz:
Tm01 = 1.086*Tz, Tp = 1.408*Tz and Tp=1.2965*Tm01
Examples
--------
>>> import wafo.spectrum.models as wsm
>>> S = wsm.Bretschneider(Hm0=6.5,Tp=10)
>>> S((0,1,2,3))
array([ 0. , 1.69350993, 0.06352698, 0.00844783])
See also
--------
Jonswap,
Torsethaugen
"""
type = 'Bretschneider'
def __init__(self, Hm0=7.0, Tp=11.0, N=5, M=4, chk_seastate=True, **kwds):
super(Bretschneider, self).__init__(Hm0, Tp)
self.N = N
self.M = M
if chk_seastate:
self.chk_seastate()
def __call__(self, wi):
""" Return Bretschnieder spectrum
"""
w = atleast_1d(wi)
if self.Hm0 > 0:
wp = 2 * pi / self.Tp
wn = w / wp
S = (self.Hm0 / 4.0) ** 2 / wp * _gengamspec(wn, self.N, self.M)
else:
S = zeros_like(w)
return S
def jonswap_peakfact(Hm0, Tp):
""" Jonswap peakedness factor, gamma, given Hm0 and Tp
Parameters
----------
Hm0 : significant wave height [m].
Tp : peak period [s]
Returns
-------
gamma : Peakedness parameter of the JONSWAP spectrum
Details
-------
A standard value for GAMMA is 3.3. However, a more correct approach is
to relate GAMMA to Hm0 and Tp:
D = 0.036-0.0056*Tp/sqrt(Hm0)
gamma = exp(3.484*(1-0.1975*D*Tp**4/(Hm0**2)))
This parameterization is based on qualitative considerations of deep water
wave data from the North Sea, see Torsethaugen et. al. (1984)
Here GAMMA is limited to 1..7.
NOTE: The size of GAMMA is the common shape of Hm0 and Tp.
Examples
--------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb
>>> Tp,Hs = plb.meshgrid(range(4,8),range(2,6))
>>> gam = wsm.jonswap_peakfact(Hs,Tp)
>>> Hm0 = plb.linspace(1,20)
>>> Tp = Hm0
>>> [T,H] = plb.meshgrid(Tp,Hm0)
>>> gam = wsm.jonswap_peakfact(H,T)
>>> v = plb.arange(0,8)
>>> Hm0 = plb.arange(1,11)
>>> Tp = plb.linspace(2,16)
>>> T,H = plb.meshgrid(Tp,Hm0)
>>> gam = wsm.jonswap_peakfact(H,T)
h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
h = plb.plot(Tp,gam.T)
h = plb.xlabel('Tp [s]')
h = plb.ylabel('Peakedness parameter')
plb.close('all')
See also
--------
jonswap
"""
Hm0, Tp = atleast_1d(Hm0, Tp)
x = Tp / sqrt(Hm0)
gam = ones_like(x)
k1 = flatnonzero(x <= 5.14285714285714)
if k1.size > 0: # limiting gamma to [1 7]
xk = x.take(k1)
D = 0.036 - 0.0056 * xk # approx 5.061*Hm0**2/Tp**4*(1-0.287*log(gam))
# gamma
gam.put(k1, minimum(exp(3.484 * (1.0 - 0.1975 * D * xk ** 4.0)), 7.0))
return gam
def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81,
output='dict'):
"""
Return Jonswap seastate from windspeed and fetch
Parameters
----------
U10 : real scalar
windspeed at 10 m above mean water surface [m/s]
fetch : real scalar
fetch [m]
method : 'hasselman73' seastate according to Hasselman et. al. 1973
'hasselman76' seastate according to Hasselman et. al. 1976
'lewis' seastate according to Lewis and Allos 1990
g : real scalar
accelaration of gravity [m/s**2]
output : 'dict' or 'list'
Returns
-------
seastate: dict where
Hm0 : significant wave height [m]
Tp : peak period [s]
gamma : jonswap peak enhancement factor.
sigmaA,
sigmaB : jonswap spectral width parameters.
Ag : jonswap alpha, normalization factor.
Example
--------
>>> import wafo.spectrum.models as wsm
>>> fetch = 10000; u10 = 10
>>> ss = wsm.jonswap_seastate(u10, fetch, output='dict')
>>> for key in sorted(ss.keys()): key, ss[key]
('Ag', 0.016257903375341734)
('Hm0', 0.51083679198275533)
('Tp', 2.7727680999585265)
('gamma', 2.4824142635861119)
('sigmaA', 0.07531733139517202)
('sigmaB', 0.09191208451225134)
>>> S = wsm.Jonswap(**ss)
>>> S.Hm0
0.51083679198275533
# Alternatively
>>> ss1 = wsm.jonswap_seastate(u10, fetch, output='list')
>>> S1 = wsm.Jonswap(*ss1)
>>> S1.Hm0
0.51083679198275533
See also
--------
Jonswap
References
----------
Lewis, A. W. and Allos, R.N. (1990)
JONSWAP's parameters: sorting out the inconscistencies.
Ocean Engng, Vol 17, No 4, pp 409-415
Hasselmann et al. (1973)
Measurements of Wind-Wave Growth and Swell Decay during the Joint
North Sea Project (JONSWAP).
Ergansungsheft, Reihe A(8), Nr. 12, Deutschen Hydrografischen Zeitschrift.
Hasselmann et al. (1976)
A parametric wave prediction model.
J. phys. oceanogr. Vol 6, pp 200-228
"""
# The following formulas are from Lewis and Allos 1990:
zeta = g * fetch / (u10 ** 2) # dimensionless fetch, Table 1
# zeta = min(zeta, 2.414655013429281e+004)
if method.startswith('h'):
if method[-1] == '3': # Hasselman et.al (1973)
A = 0.076 * zeta ** (-0.22)
# dimensionless peakfrequency, Table 1
ny = 3.5 * zeta ** (-0.33)
# dimensionless surface variance, Table 1
epsilon1 = 9.91e-8 * zeta ** 1.1
else: # Hasselman et.al (1976)
A = 0.0662 * zeta ** (-0.2)
ny = 2.84 * zeta ** (-0.3) # dimensionless peakfrequency, Table 1
# dimensionless surface variance, Eq.4
epsilon1 = 1.6e-7 * zeta
sa = 0.07
sb = 0.09
gam = 3.3
else:
A = 0.074 * zeta ** (-0.22) # Eq. 10
ny = 3.57 * zeta ** (-0.33) # dimensionless peakfrequency, Eq. 11
# dimensionless surface variance, Eq.12
epsilon1 = 3.512e-4 * A * ny ** (-4.) * zeta ** (-0.1)
sa = 0.05468 * ny ** (-0.32) # Eq. 13
sb = 0.078314 * ny ** (-0.16) # Eq. 14
gam = maximum(17.54 * zeta ** (-0.28384), 1) # Eq. 15
Tp = u10 / (ny * g) # Table 1
Hm0 = 4 * sqrt(epsilon1) * u10 ** 2. / g # Table 1
if output[0] == 'l':
return Hm0, Tp, gam, sa, sb, A
else:
return dict(Hm0=Hm0, Tp=Tp, gamma=gam, sigmaA=sa, sigmaB=sb, Ag=A)
class Jonswap(ModelSpectrum):
"""
Jonswap spectral density model
Member variables
----------------
Hm0 : significant wave height (default 7 (m))
Tp : peak period (default 11 (sec))
gamma : peakedness factor determines the concentraton
of the spectrum on the peak frequency.
Usually in the range 1 <= gamma <= 7.
default depending on Hm0, Tp, see jonswap_peakedness)
sigmaA : spectral width parameter for w<wp (default 0.07)
sigmaB : spectral width parameter for w<wp (default 0.09)
Ag : normalization factor used when gamma>1:
N : scalar defining decay of high frequency part. (default 5)
M : scalar defining spectral width around the peak. (default 4)
method : String defining method used to estimate Ag when gamma>1
'integration': Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default)
'parametric' : Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
'custom' : Ag = Ag
wnc : wc/wp normalized cut off frequency used when calculating Ag
by integration (default 6)
Parameters
----------
w : array-like
angular frequencies [rad/s]
Description
-----------
The JONSWAP spectrum is defined as
S(w) = A * Gf * G0 * wn**(-N)*exp(-N/(M*wn**M))
where
G0 = Normalizing factor related to Bretschneider form
A = Ag * (Hm0/4)**2 / wp (Normalization factor)
Gf = j**exp(-.5*((wn-1)/s)**2) (Peak enhancement factor)
wn = w/wp
wp = angular peak frequency
s = sigmaA for wn <= 1
sigmaB for 1 < wn
j = gamma, (j=1, => Bretschneider spectrum)
The JONSWAP spectrum is assumed to be especially suitable for the North
Sea, and does not represent a fully developed sea. It is a reasonable model
for wind generated sea when the seastate is in the so called JONSWAP range,
i.e., 3.6*sqrt(Hm0) < Tp < 5*sqrt(Hm0)
The relation between the peak period and mean zero-upcrossing period
may be approximated by
Tz = Tp/(1.30301-0.01698*gamma+0.12102/gamma)
Examples
---------
>>> import pylab as plb
>>> import wafo.spectrum.models as wsm
>>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
>>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
>>> w = plb.linspace(0,5)
>>> all(np.abs(S(w)-S2(w))<1.e-7)
True
h = plb.plot(w,S(w))
plb.close('all')
See also
--------
Bretschneider
Tmaspec
Torsethaugen
References
-----------
Torsethaugen et al. (1984)
Characteristica for extreme Sea States on the Norwegian continental shelf.
Report No. STF60 A84123. Norwegian Hydrodyn. Lab., Trondheim
Hasselmann et al. (1973)
Measurements of Wind-Wave Growth and Swell Decay during the Joint
North Sea Project (JONSWAP).
Ergansungsheft, Reihe A(8), Nr. 12, Deutschen Hydrografischen Zeitschrift.
"""
type = 'Jonswap'
def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09,
Ag=None, N=5, M=4, method='integration', wnc=6.0,
chk_seastate=True):
super(Jonswap, self).__init__(Hm0, Tp)
self.N = N
self.M = M
self.sigmaA = sigmaA
self.sigmaB = sigmaB
self.gamma = gamma
self.Ag = Ag
self.method = method
self.wnc = wnc
if self.gamma is None or not isfinite(self.gamma) or self.gamma < 1:
self.gamma = jonswap_peakfact(Hm0, Tp)
self._pre_calculate_ag()
if chk_seastate:
self.chk_seastate()
def _chk_extra_param(self):
Tp = self.Tp
Hm0 = self.Hm0
gam = self.gamma
outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0)
if outsideJonswapRange:
txt0 = """
Hm0=%g,Tp=%g is outside the JONSWAP range.
The validity of the spectral density is questionable.
""" % (Hm0, Tp)
warnings.warn(txt0)
if gam < 1 or 7 < gam:
txt = """
The peakedness factor, gamma, is possibly too large.
The validity of the spectral density is questionable.
"""
warnings.warn(txt)
def _localspec(self, wn):
Gf = self.peak_e_factor(wn)
return Gf * _gengamspec(wn, self.N, self.M)
def _check_parametric_ag(self, N, M, gammai):
parameters_ok = 3 <= N <= 50 or 2 <= M <= 9.5 and 1 <= gammai <= 20
if not parameters_ok:
raise ValueError('Not knowing the normalization because N, ' +
'M or peakedness parameter is out of bounds!')
if self.sigmaA != 0.07 or self.sigmaB != 0.09:
warnings.warn('Use integration to calculate Ag when ' + 'sigmaA!=0.07 or sigmaB!=0.09')
def _parametric_ag(self):
"""
Original normalization
NOTE: that Hm0**2/16 generally is not equal to intS(w)dw
with this definition of Ag if sa or sb are changed from the
default values
"""
self.method = 'parametric'
N = self.N
M = self.M
gammai = self.gamma
f1NM = 4.1 * (N - 2 * M ** 0.28 + 5.3) ** (-1.45 * M ** 0.1 + 0.96)
f2NM = ((2.2 * M ** (-3.3) + 0.57) * N ** (-0.58 * M ** 0.37 + 0.53) -
1.04 * M ** (-1.9) + 0.94)
self.Ag = (1 + f1NM * log(gammai) ** f2NM) / gammai
# if N == 5 && M == 4,
# options.Ag = (1+1.0*log(gammai).**1.16)/gammai
# options.Ag = (1-0.287*log(gammai))
# options.normalizeMethod = 'Three'
# elseif N == 4 && M == 4,
# options.Ag = (1+1.1*log(gammai).**1.19)/gammai
self._check_parametric_ag(N, M, gammai)
def _custom_ag(self):
self.method = 'custom'
if self.Ag <= 0:
raise ValueError('Ag must be larger than 0!')
def _integrate_ag(self):
# normalizing by integration
self.method = 'integration'
if self.wnc < 1.0:
raise ValueError('Normalized cutoff frequency, wnc, ' +
'must be larger than one!')
area1, unused_err1 = integrate.quad(self._localspec, 0, 1)
area2, unused_err2 = integrate.quad(self._localspec, 1, self.wnc)
area = area1 + area2
self.Ag = 1.0 / area
def _pre_calculate_ag(self):
""" PRECALCULATEAG Precalculate normalization.
"""
if self.gamma == 1:
self.Ag = 1.0
self.method = 'parametric'
elif self.Ag is not None:
self._custom_ag()
else:
norm_ag = dict(i=self._integrate_ag,
p=self._parametric_ag,
c=self._custom_ag)[self.method[0]]
norm_ag()
def peak_e_factor(self, wn):
""" PEAKENHANCEMENTFACTOR
"""
w = maximum(atleast_1d(wn), 0.0)
sab = where(w > 1, self.sigmaB, self.sigmaA)
wnm12 = 0.5 * ((w - 1.0) / sab) ** 2.0
Gf = self.gamma ** (exp(-wnm12))
return Gf
def __call__(self, wi):
""" JONSWAP spectral density
"""
w = atleast_1d(wi)
if (self.Hm0 > 0.0):
N = self.N
M = self.M
wp = 2 * pi / self.Tp
wn = w / wp
Ag = self.Ag
Hm0 = self.Hm0
Gf = self.peak_e_factor(wn)
S = ((Hm0 / 4.0) ** 2 / wp * Ag) * Gf * _gengamspec(wn, N, M)
else:
S = zeros_like(w)
return S
def phi1(wi, h, g=9.81):
""" Factor transforming spectra to finite water depth spectra.
Input
-----
w : arraylike
angular frequency [rad/s]
h : scalar
water depth [m]
g : scalar
acceleration of gravity [m/s**2]
Returns
-------
tr : arraylike
transformation factors
Example:
-------
Transform a JONSWAP spectrum to a spectrum for waterdepth = 30 m
>>> import wafo.spectrum.models as wsm
>>> S = wsm.Jonswap()
>>> w = np.arange(3.0)
>>> S(w)*wsm.phi1(w,30.0)
array([ 0. , 1.0358056 , 0.03796281])
Reference
---------
Buows, E., Gunther, H., Rosenthal, W. and Vincent, C.L. (1985)
'Similarity of the wind wave spectrum in finite depth water:
1 spectral form.'
J. Geophys. Res., Vol 90, No. C1, pp 975-986
"""
w = atleast_1d(wi)
if h == inf: # % special case infinite water depth
return ones_like(w)
k1 = w2k(w, 0, inf, g=g)[0]
dw1 = 2.0 * w / g # % dw/dk|h=inf
k2 = w2k(w, 0, h, g=g)[0]
k2h = k2 * h
den = where(k1 == 0, 1, (tanh(k2h) + k2h / cosh(k2h) ** 2.0))
dw2 = where(k1 == 0, 0, dw1 / den) # dw/dk|h=h0
return where(k1 == 0, 0, (k1 / k2) ** 3.0 * dw2 / dw1)
class Tmaspec(Jonswap):
""" JONSWAP spectrum for finite water depth
Member variables
----------------
h = water depth (default 42 [m])
g : acceleration of gravity [m/s**2]
Hm0 = significant wave height (default 7 [m])
Tp = peak period (default 11 (sec))
gamma = peakedness factor determines the concentraton
of the spectrum on the peak frequency.
Usually in the range 1 <= gamma <= 7.
default depending on Hm0, Tp, see getjonswappeakedness)
sigmaA = spectral width parameter for w<wp (default 0.07)
sigmaB = spectral width parameter for w<wp (default 0.09)
Ag = normalization factor used when gamma>1:
N = scalar defining decay of high frequency part. (default 5)
M = scalar defining spectral width around the peak. (default 4)
method = String defining method used to estimate Ag when gamma>1
'integrate' : Ag = 1/gaussq(Gf.*ggamspec(wn,N,M),0,wnc) (default)
'parametric': Ag = (1+f1(N,M)*log(gamma)^f2(N,M))/gamma
'custom' : Ag = Ag
wnc = wc/wp normalized cut off frequency used when calculating Ag
by integration (default 6)
Parameters
----------
w : array-like
angular frequencies [rad/s]
Description
------------
The evaluated spectrum is
S(w) = Sj(w)*phi(w,h)
where
Sj = jonswap spectrum
phi = modification due to water depth
The concept is based on a similarity law, and its validity is verified
through analysis of 3 data sets from: TEXEL, MARSEN projects (North
Sea) and ARSLOE project (Duck, North Carolina, USA). The data include
observations at water depths ranging from 6 m to 42 m.
Example
--------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb
>>> w = plb.linspace(0,2.5)
>>> S = wsm.Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
o=plb.plot(w,S(w))
o=plb.plot(w,S(w,h=21))
o=plb.plot(w,S(w,h=42))
plb.show()
plb.close('all')
See also
---------
Bretschneider,
Jonswap,
phi1,
Torsethaugen
References
----------
Buows, E., Gunther, H., Rosenthal, W., and Vincent, C.L. (1985)
'Similarity of the wind wave spectrum in finite depth water:
1 spectral form.'
J. Geophys. Res., Vol 90, No. C1, pp 975-986
Hasselman et al. (1973)
Measurements of Wind-Wave Growth and Swell Decay during the Joint
North Sea Project (JONSWAP).
Ergansungsheft, Reihe A(8), Nr. 12, deutschen Hydrografischen
Zeitschrift.
"""
def __init__(self, Hm0=7.0, Tp=11.0, gamma=None, sigmaA=0.07, sigmaB=0.09,
Ag=None, N=5, M=4, method='integration', wnc=6.0,
chk_seastate=True, h=42, g=9.81):
self.g = g
self.h = h
super(Tmaspec, self).__init__(Hm0, Tp, gamma, sigmaA, sigmaB, Ag, N,
M, method, wnc, chk_seastate)
self.type = 'TMA'
def phi(self, w, h=None, g=None):
if h is None:
h = self.h
if g is None:
g = self.g
return phi1(w, h, g)
def __call__(self, w, h=None, g=None):
jonswap = super(Tmaspec, self).__call__(w)
return jonswap * self.phi(w, h, g)
class Torsethaugen(ModelSpectrum):
"""
Torsethaugen double peaked (swell + wind) spectrum model
Member variables
----------------
Hm0 : significant wave height (default 7 (m))
Tp : peak period (default 11 (sec))
wnc : wc/wp normalized cut off frequency used when calculating Ag
by integration (default 6)
method : String defining method used to estimate normalization factors, Ag,
in the the modified JONSWAP spectra when gamma>1
'integrate' : Ag = 1/quad(Gf.*gengamspec(wn,N,M),0,wnc)
'parametric': Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
Parameters
----------
w : array-like
angular frequencies [rad/s]
Description
-----------
The double peaked (swell + wind) Torsethaugen spectrum is
modelled as S(w) = Ss(w) + Sw(w) where Ss and Sw are modified
JONSWAP spectrums for swell and wind peak, respectively.
The energy is divided between the two peaks according
to empirical parameters, which peak that is primary depends on parameters.
The empirical parameters are found for classes of Hm0 and Tp,
originating from a dataset consisting of 20 000 spectra divided
into 146 different classes of Hm0 and Tp. (Data measured at the
Statfjord field in the North Sea in a period from 1980 to 1989.)
The range of the measured Hm0 and Tp for the dataset
are from 0.5 to 11 meters and from 3.5 to 19 sec, respectively.
Preliminary comparisons with spectra from other areas indicate that
some of the empirical parameters are dependent on geographical location.
Thus the model must be used with care for other areas than the
North Sea and sea states outside the area where measured data
are available.
Example
-------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb
>>> w = plb.linspace(0,4)
>>> S = wsm.Torsethaugen(Hm0=6, Tp=8)
h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
See also
--------
Bretschneider
Jonswap
References
----------
Torsethaugen, K. (2004)
"Simplified Double Peak Spectral Model for Ocean Waves"
In Proc. 14th ISOPE
Torsethaugen, K. (1996)
Model for a doubly peaked wave spectrum
Report No. STF22 A96204. SINTEF Civil and Environm. Engineering, Trondheim
Torsethaugen, K. (1994)
'Model for a doubly peaked spectrum. Lifetime and fatigue strength
estimation implications.'
International Workshop on Floating Structures in Coastal zone,
Hiroshima, November 1994.
Torsethaugen, K. (1993)
'A two peak wave spectral model.'
In proceedings OMAE, Glasgow
"""
type = 'Torsethaugen'
def __init__(self, Hm0=7, Tp=11, method='integration', wnc=6, gravity=9.81,
chk_seastate=True, **kwds):
super(Torsethaugen, self).__init__(Hm0, Tp)
self.method = method
self.wnc = wnc
self.gravity = gravity
self.wind = None
self.swell = None
if chk_seastate:
self.chk_seastate()
self._init_spec()
def __call__(self, w):
""" TORSETHAUGEN spectral density
"""
return self.wind(w) + self.swell(w)
def _chk_extra_param(self):
Hm0 = self.Hm0
Tp = self.Tp
if Hm0 > 11 or Hm0 > max((Tp / 3.6) ** 2, (Tp - 2) * 12 / 11):
txt0 = """Hm0 is outside the valid range.
The validity of the spectral density is questionable"""
warnings.warn(txt0)
if Tp > 20 or Tp < 3:
txt1 = """Tp is outside the valid range.
The validity of the spectral density is questionable"""
warnings.warn(txt1)
def _init_spec(self):
""" Initialize swell and wind part of Torsethaugen spectrum
"""
monitor = 0
Hm0 = self.Hm0
Tp = self.Tp
gravity1 = self.gravity # m/s**2
min = minimum # @ReservedAssignment
max = maximum # @ReservedAssignment
# The parameter values below are found comparing the
# model to average measured spectra for the Statfjord Field
# in the Northern North Sea.
Af = 6.6 # m**(-1/3)*sec
AL = 2 # sec/sqrt(m)
Au = 25 # sec
KG = 35
KG0 = 3.5
KG1 = 1 # m
r = 0.857 # 6/7
K0 = 0.5 # 1/sqrt(m)
K00 = 3.2
M0 = 4
B1 = 2 # sec
B2 = 0.7
B3 = 3.0 # m
S0 = 0.08 # m**2*s
S1 = 3 # m
# Preliminary comparisons with spectra from other areas indicate that
# the parameters on the line below can be dependent on geographical
# location
A10 = 0.7
A1 = 0.5
A20 = 0.6
A2 = 0.3
A3 = 6
Tf = Af * (Hm0) ** (1.0 / 3.0)
Tl = AL * sqrt(Hm0) # lower limit
Tu = Au # upper limit
# Non-dimensional scales
# New call pab April 2005
El = min(max((Tf - Tp) / (Tf - Tl), 0), 1) # wind sea
Eu = min(max((Tp - Tf) / (Tu - Tf), 0), 1) # Swell
if Tp < Tf: # Wind dominated seas
# Primary peak (wind dominated)
Nw = K0 * sqrt(Hm0) + K00 # high frequency exponent
Mw = M0 # spectral width exponent
Rpw = min((1 - A10) * exp(-(El / A1) ** 2) + A10, 1)
Hpw = Rpw * Hm0 # significant waveheight wind
Tpw = Tp # primary peak period
# peak enhancement factor
gammaw = KG * (1 + KG0 * exp(-Hm0 / KG1)) * \
(2 * pi / gravity1 * Rpw * Hm0 / (Tp ** 2)) ** r
gammaw = max(gammaw, 1)
# Secondary peak (swell)
Ns = Nw # high frequency exponent
Ms = Mw # spectral width exponent
Rps = sqrt(1.0 - Rpw ** 2.0)
Hps = Rps * Hm0 # significant waveheight swell
Tps = Tf + B1
gammas = 1.0
if monitor:
if Rps > 0.1:
print(' Spectrum for Wind dominated sea')
else:
print(' Spectrum for pure wind sea')
else: # swell dominated seas
# Primary peak (swell)
Ns = K0 * sqrt(Hm0) + K00 # high frequency exponent
Ms = M0 # spectral width exponent
Rps = min((1 - A20) * exp(-(Eu / A2) ** 2) + A20, 1)
Hps = Rps * Hm0 # significant waveheight swell
Tps = Tp # primary peak period
# peak enhancement factor
gammas = KG * (1 + KG0 * exp(-Hm0 / KG1)) * \
(2 * pi / gravity1 * Hm0 / (Tf ** 2)) ** r * (1 + A3 * Eu)
gammas = max(gammas, 1)
# Secondary peak (wind)
Nw = Ns # high frequency exponent
Mw = M0 * (1 - B2 * exp(-Hm0 / B3)) # spectral width exponent
Rpw = sqrt(1 - Rps ** 2)
Hpw = Rpw * Hm0 # significant waveheight wind
C = (Nw - 1) / Mw
B = Nw / Mw
G0w = B ** C * Mw / sp.gamma(C) # normalizing factor
# G0w = exp(C*log(B)+log(Mw)-gammaln(C))
# G0w = Mw/((B)**(-C)*gamma(C))
if Hpw > 0:
Tpw = (16 * S0 * (1 - exp(-Hm0 / S1)) * (0.4) **
Nw / (G0w * Hpw ** 2)) ** (-1.0 / (Nw - 1.0))
else:
Tpw = inf
# Tpw = max(Tpw,2.5)
gammaw = 1
if monitor:
if Rpw > 0.1:
print(' Spectrum for swell dominated sea')
else:
print(' Spectrum for pure swell sea')
if monitor:
if (3.6 * sqrt(Hm0) <= Tp & Tp <= 5 * sqrt(Hm0)):
print(' Jonswap range')
print('Hm0 = %g' % Hm0)
print('Ns, Ms = %g, %g Nw, Mw = %g, %g' % (Ns, Ms, Nw, Mw))
print('gammas = %g gammaw = %g' % (gammas, gammaw))
print('Rps = %g Rpw = %g' % (Rps, Rpw))
print('Hps = %g Hpw = %g' % (Hps, Hpw))
print('Tps = %g Tpw = %g' % (Tps, Tpw))
# G0s=Ms/((Ns/Ms)**(-(Ns-1)/Ms)*gamma((Ns-1)/Ms )) #normalizing factor
self.wind = Jonswap(Hm0=Hpw, Tp=Tpw, gamma=gammaw, N=Nw, M=Mw,
method=self.method, chk_seastate=False)
self.swell = Jonswap(Hm0=Hps, Tp=Tps, gamma=gammas, N=Ns, M=Ms,
method=self.method, chk_seastate=False)
class McCormick(Bretschneider):
""" McCormick spectral density model
Member variables
----------------
Hm0 = significant wave height (default 7 (m))
Tp = peak period (default 11 (sec))
Tz = zero-down crossing period (default 0.8143*Tp)
M = scalar defining spectral width around the peak.
(default depending on Tp and Tz)
Parameters
----------
w : array-like
angular frequencies [rad/s]
Description
-----------
The McCormick spectrum parameterization is a modification of the
Bretschneider spectrum and defined as
S(w) = (M+1)*(Hm0/4)^2/wp*(wp./w)^(M+1)*exp(-(M+1)/M*(wp/w)^M)
where
Tp/Tz=(1+1/M)^(1/M)/gamma(1+1/M)
Example:
--------
>>> import wafo.spectrum.models as wsm
>>> S = wsm.McCormick(Hm0=6.5,Tp=10)
>>> S(range(4))
array([ 0. , 1.87865908, 0.15050447, 0.02994663])
See also
--------
Bretschneider
Jonswap,
Torsethaugen
References:
-----------
M.E. McCormick (1999)
"Application of the Generic Spectral Formula to Fetch-Limited Seas"
Marine Technology Society, Vol 33, No. 3, pp 27-32
"""
type = 'McCormick'
def __init__(self, Hm0=7, Tp=11, Tz=None, M=None, chk_seastate=True):
if Tz is None:
Tz = 0.8143 * Tp
self.Tz = Tz
if M is None and Hm0 > 0:
self._TpdTz = Tp / Tz
M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5)
N = M + 1.0
super(McCormick, self).__init__(Hm0, Tp, N, M, chk_seastate)
def _localoptfun(self, x):
# LOCALOPTFUN Local function to optimize.
y = 1.0 + x
return (y ** (x) / sp.gamma(y) - self._TpdTz) ** 2.0
class OchiHubble(ModelSpectrum):
""" OchiHubble bimodal spectral density model.
Member variables
----------------
Hm0 : significant wave height (default 7 (m))
par : integer defining the parametrization (default 0)
0 : The most probable spectrum
1,2,...10 : gives 95% Confidence spectra
The OchiHubble bimodal spectrum is modelled as
S(w) = Ss(w) + Sw(w) where Ss and Sw are modified Bretschneider
spectra for swell and wind peak, respectively.
The OH spectrum is a six parameter spectrum, all functions of Hm0.
The values of these parameters are determined from a analysis of data
obtained in the North Atlantic. The source of the data is the same as
that for the development of the Pierson-Moskowitz spectrum, but
analysis is carried out on over 800 spectra including those in
partially developed seas and those having a bimodal shape. From a
statistical analysis of the data, a family of wave spectra consisting
of 11 members is generated for a desired sea severity (Hm0) with the
coefficient of 0.95.
A significant advantage of using a family of spectra for design of
marine systems is that one of the family members yields the largest
response such as motions or wave induced forces for a specified sea
severity, while another yields the smallest response with confidence
coefficient of 0.95.
Examples
--------
>>> import wafo.spectrum.models as wsm
>>> S = wsm.OchiHubble(par=2)
>>> S(range(4))
array([ 0. , 0.90155636, 0.04185445, 0.00583207])
See also
--------
Bretschneider,
Jonswap,
Torsethaugen
References:
----------
Ochi, M.K. and Hubble, E.N. (1976)
'On six-parameter wave spectra.'
In Proc. 15th Conf. Coastal Engng., Vol.1, pp301-328
"""
type = 'Ochi Hubble'
def __init__(self, Hm0=7, par=1, chk_seastate=True):
super(OchiHubble, self).__init__(Hm0, Tp=1)
self.par = par
self.wind = None
self.swell = None
if chk_seastate:
self.chk_seastate()
self._init_spec()
def __call__(self, w):
return self.wind(w) + self.swell(w)
def _init_spec(self):
hp = array([[0.84, 0.54],
[0.84, 0.54],
[0.84, 0.54],
[0.84, 0.54],
[0.84, 0.54],
[0.95, 0.31],
[0.65, 0.76],
[0.90, 0.44],
[0.77, 0.64],
[0.73, 0.68],
[0.92, 0.39]])
wa = array([[0.7, 1.15],
[0.93, 1.5],
[0.41, 0.88],
[0.74, 1.3],
[0.62, 1.03],
[0.70, 1.50],
[0.61, 0.94],
[0.81, 1.60],
[0.54, 0.61],
[0.70, 0.99],
[0.70, 1.37]])
wb = array([[0.046, 0.039],
[0.056, 0.046],
[0.016, 0.026],
[0.052, 0.039],
[0.039, 0.030],
[0.046, 0.046],
[0.039, 0.036],
[0.052, 0.033],
[0.039, 0.000],
[0.046, 0.039],
[0.046, 0.039]])
Lpar = array([[3.00, 1.54, -0.062],
[3.00, 2.77, -0.112],
[2.55, 1.82, -0.089],
[2.65, 3.90, -0.085],
[2.60, 0.53, -0.069],
[1.35, 2.48, -0.102],
[4.95, 2.48, -0.102],
[1.80, 2.95, -0.105],
[4.50, 1.95, -0.082],
[6.40, 1.78, -0.069],
[0.70, 1.78, -0.069]])
Hm0 = self.Hm0
Lpari = Lpar[self.par]
Li = hstack((Lpari[0], Lpari[1] * exp(Lpari[2] * Hm0)))
Hm0i = hp[self.par] * Hm0
Tpi = 2 * pi * exp(wb[self.par] * Hm0) / wa[self.par]
Ni = 4 * Li + 1
Mi = [4, 4]
self.swell = Bretschneider(Hm0=Hm0i[0], Tp=Tpi[0], N=Ni[0], M=Mi[0])
self.wind = Bretschneider(Hm0=Hm0i[1], Tp=Tpi[1], N=Ni[1], M=Mi[1])
def _chk_extra_param(self):
if self.par < 0 or 10 < self.par:
raise ValueError('Par must be an integer from 0 to 10!')
class Wallop(Bretschneider):
"""Wallop spectral density model.
Member variables
----------------
Hm0 = significant wave height (default 7 (m))
Tp = peak period (default 11 (sec))
N = shape factor, i.e. slope for the high frequency
% part (default depending on Hm0 and Tp, see below)
Parameters
----------
w : array-like
angular frequencies [rad/s]
Description
-----------
The WALLOP spectrum parameterization is a modification of the Bretschneider
spectrum and defined as
S(w) = A * G0 * wn**(-N)*exp(-N/(4*wn**4))
where
G0 = Normalizing factor related to Bretschneider form
A = (Hm0/4)^2 / wp (Normalization factor)
wn = w/wp
wp = 2*pi/Tp, angular peak frequency
N = abs((log(2*pi^2)+2*log(Hm0/4)-2*log(Lp))/log(2))
Lp = wave length corresponding to the peak frequency, wp.
If N=5 it becomes the same as the JONSWAP spectrum with
peak enhancement factor gamma=1 or the Bretschneider
(Pierson-Moskowitz) spectrum.
Example:
--------
>>> import wafo.spectrum.models as wsm
>>> S = wsm.Wallop(Hm0=6.5, Tp=10)
>>> S(range(4))
array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03,
7.72996150e-05])
See also
--------
Bretschneider
Jonswap,
Torsethaugen
References:
-----------
Huang, N.E., Long, S.R., Tung, C.C, Yuen, Y. and Bilven, L.F. (1981)
"A unified two parameter wave spectral model for a generous sea state"
J. Fluid Mechanics, Vol.112, pp 203-224
"""
type = 'Wallop'
def __init__(self, Hm0=7, Tp=11, N=None, chk_seastate=True):
M = 4
if N is None:
wp = 2. * pi / Tp
kp = w2k(wp, 0, inf)[0] # wavenumber at peak frequency
Lp = 2. * pi / kp # wave length at the peak frequency
N = np.abs((log(2. * pi ** 2.) + 2 * log(Hm0 / 4) -
2.0 * log(Lp)) / log(2))
super(Wallop, self).__init__(Hm0, Tp, N, M, chk_seastate)
class Spreading(object):
"""
Directional spreading function.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
wc : real scalar
cut over frequency
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
Member variables
----------------
type : string (default 'cos-2s')
type of spreading function, see options below
'cos-2s' : N(S)*[cos((theta-theta0)/2)]**(2*S) (0 < S)
'Box-car' : N(A)*I( -A < theta-theta0 < A) (0 < A < pi)
'von-Mises' : N(K)*exp(K*cos(theta-theta0)) (0 < K)
'Poisson' : N(X)/(1-2*X*cos(theta-theta0)+X**2) (0 < X < 1)
'sech-2' : N(B)*sech(B*(theta-theta0))**2 (0 < B)
'wrapped-normal':
[1 + 2*sum exp(-(n*D1)^2/2)*cos(n*(theta-theta0))]/(2*pi) (0 < D1)
(N(.) = normalization factor)
(the first letter is enough for unique identification)
theta0 : callable, matrix or a scalar
defines average direction given in radians at every angular frequency.
(length 1 or length == length(wn)) (default 0)
method : string or integer
Defines function used for direcional spreading parameter:
0, None : S(wn) = s_a, frequency independent
1, 'mitsuyasu': S(wn) frequency dependent (default)
where S(wn) = s_a *(wn)**m_a, for wn_lo <= wn < wn_c
= s_b *(wn)**m_b, for wn_c <= wn < wn_up
= 0 otherwise
2, 'donelan' : B(wn) frequency dependent
3, 'banner' : B(wn) frequency dependent
where B(wn) = S(wn) for wn_lo <= wn < wn_up
= s_b*wn_up**m_b, for wn_up <= wn and method = 2
= sc*F(wn) for wn_up <= wn and method = 3
where F(wn) = 10^(-0.4+0.8393*exp(-0.567*log(wn^2))) and
sc is scalefactor to make the spreading funtion continous.
wn_lo, wn_c, wn_up: real scalars (default 0, 1, inf)
limits used in the function defining the directional spreading
parameter, S() or B() defined above.
wn_c is the normalized cutover frequency
s_a, s_b : real scalars
maximum spread parameters (default [15 15])
m_a, m_b : real scalars
shape parameters (default [5 -2.5])
SPREADING return a Directional spreading function.
Here the S- or B-parameter, of the COS-2S and SECH-2 spreading function,
respectively, is used as a measure of spread. All the parameters of the
other distributions are related to this parameter through the first Fourier
coefficient, R1, of the directional distribution as follows:
R1 = S/(S+1) or S = R1/(1-R1).
where
Box-car spreading : R1 = sin(A)/A
Von Mises spreading: R1 = besseli(1,K)/besseli(0,K),
Poisson spreading : R1 = X
sech-2 spreading : R1 = pi/(2*B*sinh(pi/(2*B))
Wrapped Normal : R1 = exp(-D1^2/2)
A value of S = 15 corresponds to
'box' : A=0.62, 'sech-2' : B=0.89
'von-mises' : K=8.3, 'poisson': X=0.94
'wrapped-normal': D=0.36
The COS2S is the most frequently used spreading in engineering practice.
Apart from the current meter/pressure cell data in WADIC all
instruments seem to support the 'cos2s' distribution for heavier sea
states, (Krogstad and Barstow, 1999). For medium sea states
a spreading function between COS2S and POISSON seem appropriate,
while POISSON seems appropriate for swell.
For the COS2S Mitsuyasu et al. parameterized SPa = SPb =
11.5*(U10/Cp) where Cp = g/wp is the deep water phase speed at wp and
U10 the wind speed at reference height 10m. Hasselman et al. (1980)
parameterized mb = -2.33-1.45*(U10/Cp-1.17).
Mitsuyasu et al. (1975) showed that SP for wind waves varies from
5 to 30 being a function of dimensionless wind speed.
However, Goda and Suzuki (1975) proposed SP = 10 for wind waves, SP = 25
for swell with short decay distance and SP = 75 for long decay distance.
Compared to experiments Krogstad et al. (1998) found that m_a = 5 +/- _EPS
and that -1< m_b < -3.5.
Values given in the litterature: [s_a s_b m_a m_b wn_lo wn_c wn_up]
(Mitsuyasu: s_a == s_b) (cos-2s) [15 15 5 -2.5 0 1 3 ]
(Hasselman: s_a ~= s_b) (cos-2s) [6.97 9.77 4.06 -2.3 0 1.05 3 ]
(Banner : s_a ~= s_b) (sech2) [2.61 2.28 1.3 -1.3 0.56 0.95 1.6]
Examples
--------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb
>>> D = wsm.Spreading('cos2s',s_a=10.0)
# Make directionale spectrum
>>> S = wsm.Jonswap().tospecdata()
>>> SD = D.tospecdata2d(S)
>>> w = plb.linspace(0,3,257)
>>> theta = plb.linspace(-pi,pi,129)
# Make frequency dependent direction spreading
>>> theta0 = lambda w: w*plb.pi/6.0
>>> D2 = wsm.Spreading('cos2s',theta0=theta0)
h = SD.plot()
t = plb.contour(D(theta,w)[0].squeeze())
t = plb.contour(D2(theta,w)[0])
# Plot all spreading functions
alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm')
for ix in range(len(alltypes)):
... D3 = wsm.Spreading(alltypes[ix])
... t = plb.figure(ix)
... t = plb.contour(D3(theta,w)[0])
... t = plb.title(alltypes[ix])
plb.close('all')
See also
--------
mkdspec, plotspec, spec2spec
References
---------
Krogstad, H.E. and Barstow, S.F. (1999)
"Directional Distributions in Ocean Wave Spectra"
Proceedings of the 9th ISOPE Conference, Vol III, pp. 79-86
Goda, Y. (1999)
"Numerical simulation of ocean waves for statistical analysis"
Marine Tech. Soc. Journal, Vol. 33, No. 3, pp 5--14
Banner, M.L. (1990)
"Equilibrium spectra of wind waves."
J. Phys. Ocean, Vol 20, pp 966--984
Donelan M.A., Hamilton J, Hui W.H. (1985)
"Directional spectra of wind generated waves."
Phil. Trans. Royal Soc. London, Vol A315, pp 387--407
Hasselmann D, Dunckel M, Ewing JA (1980)
"Directional spectra observed during JONSWAP."
J. Phys. Ocean, Vol.10, pp 1264--1280
Mitsuyasu, H, et al. (1975)
"Observation of the directional spectrum of ocean waves using a
coverleaf buoy."
J. Physical Oceanography, Vol.5, No.4, pp 750--760
Some of this might be included in help header:
cos-2s:
NB! The generally strong frequency dependence in directional spread
makes it questionable to run load tests of ships and structures with a
directional spread independent of frequency (Krogstad and Barstow, 1999).
"""
# Parameterization of B
# def = 2 Donelan et al freq. parametrization for 'sech2'
# def = 3 Banner freq. parametrization for 'sech2'
# (spa ~= spb) (sech-2) [2.61 2.28 1.3 -1.3 0.56 0.95 1.6]
#
def __init__(self, type='cos-2s', theta0=0, # @ReservedAssignment
method='mitsuyasu', s_a=15., s_b=15., m_a=5., m_b=-2.5,
wn_lo=0.0, wn_c=1., wn_up=inf):
self.type = type
self.theta0 = theta0
self.method = method
self.s_a = s_a
self.s_b = s_b
self.m_a = m_a
self.m_b = m_b
self.wn_lo = wn_lo
self.wn_c = wn_c
self.wn_up = wn_up
self._spreadfun = dict(c=self.cos2s, b=self.box, m=self.mises,
v=self.mises,
p=self.poisson, s=self.sech2, w=self.wrap_norm)
self._fourierdispatch = dict(b=self.fourier2a, m=self.fourier2k,
v=self.fourier2k,
p=self.fourier2x, s=self.fourier2b,
w=self.fourier2d)
@property
def method(self):
return self._method
@method.setter
def method(self, method):
methods = {'n': None, 'm': 'mitsuyasu', 'd': 'donelan', 'b': 'banner',
0: None, 1: 'mitsuyasu', 2: 'donelan', 3: 'banner',
None: None}
m = method if not isinstance(method, str) else method[0].lower()
try:
self._method = methods[m]
except KeyError:
msg = 'Unknown method. Got {}, but expected one of {}!'
raise ValueError(msg.format(method, str(methods.keys())))
def __call__(self, theta, w=1, wc=1):
spreadfun = self._spreadfun[self.type[0]]
return spreadfun(theta, w, wc)
def _normalize_angle(self, wn, theta, th0):
Nt0 = th0.size
Nw = wn.size
isFreqDepDir = (Nt0 == Nw)
if isFreqDepDir:
# frequency dependent spreading and/or
# frequency dependent direction
# make sure -pi<=TH<pi
TH = mod(theta[:, newaxis] - th0[newaxis, :] + pi, 2 * pi) - pi
elif Nt0 != 1:
raise ValueError(
'The length of theta0 must equal to 1 or the length of w')
else:
TH = mod(theta - th0 + pi, 2 * pi) - pi # make sure -pi<=TH<pi
if self.method is not None: # frequency dependent spreading
TH = TH[:, newaxis]
return TH
def _get_main_direction(self, wn):
if hasattr(self.theta0, '__call__'):
return self.theta0(wn.flatten())
return atleast_1d(self.theta0).flatten()
def chk_input(self, theta, w=1, wc=1):
""" CHK_INPUT
CALL [s_par,TH,phi0,Nt] = inputchk(theta,w,wc)
"""
wn = atleast_1d(w / wc)
theta = theta.ravel()
Nt = len(theta)
# Make sure theta is from -pi to pi
phi0 = 0.0
theta = mod(theta + pi, 2 * pi) - pi
theta0 = self._get_main_direction(wn)
TH = self._normalize_angle(wn, theta, theta0)
s = self.spread_parameter_s(wn)
return s, TH, phi0, Nt
def cos2s(self, theta, w=1, wc=1): # [D, phi0] =
""" COS2S spreading function
cos2s(theta,w) = N(S)*[cos((theta-theta0)/2)]^(2*S) (0 < S)
where N() is a normalization factor and S is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
S, TH, phi0 = self.chk_input(theta, w, wc)[:3]
gammaln = sp.gammaln
D = (exp(gammaln(S + 1) - gammaln(S + 1.0 / 2.0)) / (2 * sqrt(pi))) * \
cos(TH / 2.0) ** (2.0 * S)
return D, phi0
def poisson(self, theta, w=1, wc=1): # [D,phi0] =
""" POISSON spreading function
poisson(theta,w) = N(X)/(1-2*X*cos(theta-theta0)+X^2) (0 < X < 1)
where N() is a normalization factor and X is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
X, TH, phi0 = self.chk_input(theta, w, wc)[:3]
D = (1 - X ** 2.) / (1. - (2. * cos(TH) - X) * X) / (2. * pi)
return D, phi0
def wrap_norm(self, theta, w=1, wc=1):
""" Wrapped Normal spreading function
wnormal(theta,w) = N(D1)*[1 +
2*sum exp(-(n*D1)^2/2)*cos(n*(theta-theta0))] (0 < D1)
where N() is a normalization factor and D1 is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
par, TH, phi0, Nt = self.chk_input(theta, w, wc)
D1 = par ** 2. / 2.
ix = arange(1, Nt)
ix2 = ix ** 2
Nd2 = D1.size
Fcof = vstack((ones((1, Nd2)) / 2, exp(-ix2[:, newaxis] * D1))) / pi
cor = exp(1j * ix[:, newaxis] * TH[0, :])
# correction term to get
Pcor = vstack((ones((1, TH.shape[1])), cor))
# the correct integration limits
Fcof = Fcof * Pcor.conj()
D = real(fft(Fcof, axis=0))
D[D < 0] = 0
return D, phi0
def sech2(self, theta, w=1, wc=1):
"""SECH2 directonal spreading function
sech2(theta,w) = N(B)*0.5*B*sech(B*(theta-theta0))^2 (0 < B)
where N() is a normalization factor and X is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
B, TH, phi0 = self.chk_input(theta, w, wc)[:3]
NB = tanh(pi * B) # % Normalization factor.
NB = where(NB == 0, 1.0, NB) # Avoid division by zero
D = 0.5 * B * sech(B * TH) ** 2. / NB
return D, phi0
def mises(self, theta, w=1, wc=1):
"""Mises spreading function
mises(theta,w) = N(K)*exp(K*cos(theta-theta0)) (0 < K)
where N() is a normalization factor and K is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : arrays
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
K, TH, phi0 = self.chk_input(theta, w, wc)[:3]
D = exp(K * (cos(TH) - 1.)) / (2 * pi * sp.ive(0, K))
return D, phi0
def box(self, theta, w=1, wc=1):
""" Box car spreading function
box(theta,w) = N(A)*I( -A < theta-theta0 < A) (0 < A < pi)
where N() is a normalization factor and A is the spreading parameter
possibly dependent on w.
Parameters
----------
theta, w : vectors
angles and angular frequencies given in radians and rad/s,
respectively. Lenghts are Nt and Nw.
Returns
-------
D : 2D array
Directonal spreading function. size Nt X Nw.
The principal direction of D is always along the x-axis.
phi0 : real scalar
Parameter defining the actual principal direction of D.
"""
A, TH, phi0 = self.chk_input(theta, w, wc)[:3]
D = ((-A <= TH) & (TH <= A)) / (2. * A)
return D, phi0
# Local sub functions
def fourier2distpar(self, r1):
""" Fourier coefficients to distribution parameter
Parameters
----------
r1 = corresponding fourier coefficient.
type = string defining spreading function
'box'
'mises'
'poisson'
'sech2'
'wnormal'
Returns
x = distribution parameter
The S-parameter of the COS-2S spreading function is used as a measure
of spread in MKSPREADING. All the parameters of the other
distributions are related to this S-parameter through the first
Fourier coefficient, R1, of the directional distribution as follows:
R1 = S/(S+1) or S = R1/(1-R1).
where
Box-car spreading : R1 = sin(A)/A
Von Mises spreading: R1 = besseli(1,K)/besseli(0,K),
Poisson spreading : R1 = X
sech-2 spreading : R1 = pi/(2*B*sinh(pi/(2*B))
Wrapped Normal : R1 = exp(-D1^2/2)
"""
fourierfun = self._fourierdispatch.get(self.type[0])
return fourierfun(r1)
@staticmethod
def fourier2x(r1):
""" Returns the solution of r1 = x.
"""
X = r1
if np.any(X >= 1):
raise ValueError('POISSON spreading: X value must be less than 1')
return X
@staticmethod
def fourier2a(r1):
""" Returns the solution of R1 = sin(A)/A.
"""
A0 = flipud(linspace(0, pi + 0.1, 1025))
funA = interp1d(sinc(A0 / pi), A0)
A0 = funA(r1.ravel())
A = asarray(A0)
# Newton-Raphson
da = ones_like(r1)
max_count = 100
ix = flatnonzero(A)
for unused_iy in range(max_count):
Ai = A[ix]
da[ix] = (sin(Ai) - Ai * r1[ix]) / (cos(Ai) - r1[ix])
Ai = Ai - da[ix]
# Make sure that the current guess is larger than zero.
A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0)
ix = flatnonzero(
(np.abs(da) > sqrt(_EPS) * np.abs(A)) * (np.abs(da) > sqrt(_EPS)))
if ix.size == 0:
if np.any(A > pi):
raise ValueError(
'BOX-CAR spreading: The A value must be less than pi')
return A.clip(min=1e-16, max=pi)
warnings.warn('Newton raphson method did not converge.')
return A.clip(min=1e-16) # Avoid division by zero
@staticmethod
def fourier2k(r1):
"""
Returns the solution of R1 = besseli(1,K)/besseli(0,K),
"""
def fun0(x):
return sp.ive(1, x) / sp.ive(0, x)
K0 = hstack((linspace(0, 10, 513), linspace(10.00001, 100)))
funK = interp1d(fun0(K0), K0)
K0 = funK(r1.ravel())
k1 = flatnonzero(isnan(K0))
if (k1.size > 0):
K0[k1] = 0.0
K0[k1] = K0.max()
ix0 = flatnonzero(r1 != 0.0)
K = zeros_like(r1)
for ix in ix0:
K[ix] = optimize.fsolve(lambda x: fun0(x) - r1[ix], K0[ix])
return K
def fourier2b(self, r1):
""" Returns the solution of R1 = pi/(2*B*sinh(pi/(2*B)).
"""
B0 = hstack((linspace(_EPS, 5, 513), linspace(5.0001, 100)))
funB = interp1d(self._r1ofsech2(B0), B0)
B0 = funB(r1.ravel())
k1 = flatnonzero(isnan(B0))
if (k1.size > 0):
B0[k1] = 0.0
B0[k1] = max(B0)
ix0 = flatnonzero(r1 != 0.0)
B = zeros_like(r1)
def fun(x):
return 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix]
for ix in ix0:
B[ix] = np.abs(optimize.fsolve(fun, B0[ix]))
return B
def fourier2d(self, r1):
""" Returns the solution of R1 = exp(-D**2/2).
"""
r = clip(r1, 0., 1.0)
return where(r <= 0, inf, sqrt(-2.0 * log(r)))
def _init_frequency_dependent_spreading(self, wn):
wn_lo, wn_up = self.wn_lo, self.wn_up
wn_c = self.wn_c
spa, spb = self.s_a, self.s_b
ma, mb = self.m_a, self.m_b
# Mitsuyasu et. al and Hasselman et. al parametrization of
# frequency dependent spreading
s = where(wn <= wn_c, spa * wn ** ma, spb * wn ** mb)
s[wn <= wn_lo] = 0.0
return s, spb, wn_up, mb
def _donelan_spread(self, wn):
# Donelan et. al. parametrization for B in SECH-2
s, spb, wn_up, mb = self._init_frequency_dependent_spreading(wn)
k = flatnonzero(wn_up < wn)
s[k] = spb * (wn_up) ** mb
# Convert to S-paramater in COS-2S distribution
r1 = self.r1ofsech2(s)
s = r1 / (1. - r1)
return s
def _banner_spread(self, wn):
# Donelan et. al. parametrization for B in SECH-2
s, spb, wn_up, mb = self._init_frequency_dependent_spreading(wn)
k = flatnonzero(wn_up < wn)
# Banner parametrization for B in SECH-2
s3m = spb * (wn_up) ** mb
s3p = self._donelan(wn_up)
# Scale so that parametrization will be continous
scale = s3m / s3p
s[k] = scale * self.donelan(wn[k])
r1 = self.r1ofsech2(s)
# Convert to S-paramater in COS-2S distribution
s = r1 / (1. - r1)
return s
def _mitsuyasu_spread(self, wn):
s, _spb, wn_up, _mb = self._init_frequency_dependent_spreading(wn)
k = flatnonzero(wn_up < wn)
s[k] = 0
return s
def _frequency_independent_spread(self, _wn):
"""
no frequency dependent spreading,
but possible frequency dependent direction
"""
return atleast_1d(self.s_a)
def spread_parameter_s(self, wn):
""" Return spread parameter, S, equivalent for the COS2S function
Parameters
----------
wn : array_like
normalized frequencies.
Returns
-------
S : ndarray
spread parameter of COS2S functions
"""
spread = dict(b=self._banner_spread,
d=self._donelan_spread,
m=self._mitsuyasu_spread
).get(self.method[0],
self._frequency_independent_spread)
s = spread(wn)
if np.any(s < 0):
raise ValueError('The COS2S spread parameter, S(w), ' +
'value must be larger than 0')
if self.type[0] == 'c': # cos2s
s_par = s
else:
# First Fourier coefficient of the directional spreading function.
r1 = np.abs(s / (s + 1))
# Find distribution parameter from first Fourier coefficient.
s_par = self.fourier2distpar(r1)
if self.method is not None:
s_par = s_par[newaxis, :]
return s_par
@staticmethod
def _donelan(wn):
""" High frequency decay of B of sech2 paramater
"""
return 10.0 ** (-0.4 + 0.8393 * exp(-0.567 * log(wn ** 2)))
@staticmethod
def _r1ofsech2(B):
""" R1OFSECH2 Computes R1 = pi./(2*B.*sinh(pi./(2*B)))
"""
realmax = finfo(float).max
tiny = 1. / realmax
x = clip(2. * B, tiny, realmax)
xk = pi / x
return where(x < 100., xk / sinh(xk),
-2. * xk / (exp(xk) * expm1(-2. * xk)))
@staticmethod
def _check_theta(theta):
L = np.abs(theta[-1] - theta[0])
if np.abs(L - 2 * np.pi) > _EPS:
raise ValueError('theta must cover all angles -pi -> pi')
nt = len(theta)
if nt < 40:
warnings.warn('Number of angles is less than 40. ' +
'Spreading too sparsely sampled!')
def tospecdata2d(self, specdata, theta=None, wc=0.52, nt=51):
"""
MKDSPEC Make a directional spectrum
frequency spectrum times spreading function
CALL: Snew=mkdspec(S,D,plotflag)
Snew = directional spectrum (spectrum struct)
S = frequency spectrum (spectrum struct)
(default jonswap)
D = spreading function (special struct)
(default spreading([],'cos2s'))
Creates a directional spectrum through multiplication of a frequency
spectrum and a spreading function: S(w,theta)=S(w)*D(w,theta)
The spreading structure must contain the following fields:
.S (size [np 1] or [np nf]) and .theta (length np)
optional fields: .w (length nf), .note (memo) .phi (rotation-azymuth)
NB! S.w and D.w (if any) must be identical.
Example
-------
>>> import wafo.spectrum.models as wsm
>>> S = wsm.Jonswap().tospecdata()
>>> D = wsm.Spreading('cos2s')
>>> SD = D.tospecdata2d(S)
h = SD.plot()
See also spreading, rotspec, jonswap, torsethaugen
"""
if theta is None:
theta = np.linspace(-np.pi, np.pi, nt)
self._check_theta(theta)
w = specdata.args
S = specdata.data
D, phi0 = self(theta, w=w, wc=wc)
if D.ndim != 2: # frequency dependent spreading
D = D[:, None]
SD = D * S[None, :]
Snew = SpecData2D(SD, (w, theta), type='dir',
freqtype=specdata.freqtype)
Snew.tr = specdata.tr
Snew.h = specdata.h
Snew.phi = phi0
Snew.norm = specdata.norm
# Snew.note = specdata.note + ', spreading: %s' % self.type
return Snew
def _test_some_spectra():
S = Jonswap()
w = arange(3.0)
S(w) * phi1(w, 30.0)
S1 = S.tospecdata(w)
S1.plot()
import pylab as plb
w = plb.linspace(0, 2.5)
S = Tmaspec(h=10, gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
plb.plot(w, S(w))
plb.plot(w, S(w, h=21))
plb.plot(w, S(w, h=42))
plb.show()
plb.close('all')
w, th = plb.ogrid[0:4, 0:6]
k, k2 = w2k(w, th)
plb.plot(w, k, w, k2)
plb.show()
plb.close('all')
w = plb.linspace(0, 2, 100)
S = Torsethaugen(Hm0=6, Tp=8)
plb.plot(w, S(w), w, S.wind(w), w, S.swell(w))
S1 = Jonswap(Hm0=7, Tp=11, gamma=1)
w = plb.linspace(0, 2, 100)
plb.plot(w, S1(w))
plb.show()
plb.close('all')
Hm0 = plb.arange(1, 11)
Tp = plb.linspace(2, 16)
T, H = plb.meshgrid(Tp, Hm0)
gam = jonswap_peakfact(H, T)
plb.plot(Tp, gam.T)
plb.xlabel('Tp [s]')
plb.ylabel('Peakedness parameter')
Hm0 = plb.linspace(1, 20)
Tp = Hm0
[T, H] = plb.meshgrid(Tp, Hm0)
gam = jonswap_peakfact(H, T)
v = plb.arange(0, 8)
plb.contourf(Tp, Hm0, gam, v)
plb.colorbar()
plb.show()
plb.close('all')
def _test_spreading():
import pylab as plb
pi = plb.pi
w = plb.linspace(0, 3, 257)
theta = plb.linspace(-pi, pi, 129)
D2 = Spreading('cos2s', theta0=lambda w: w * plb.pi / 6.0)
d1 = D2(theta, w)[0]
plb.contour(d1.squeeze())
pi = plb.pi
D = Spreading('wrap_norm', s_a=10.0)
w = plb.linspace(0, 3, 257)
theta = plb.linspace(-pi, pi, 129)
d1 = D(theta, w)
plb.contour(d1[0])
plb.show()
def test_docstrings():
import doctest
print('Testing docstrings in %s' % __file__)
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
def main():
if False: # True: #
_test_some_spectra()
else:
test_docstrings()
if __name__ == '__main__':
main()
Reference in New Issue View Git Blame Copy Permalink