danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

updates

Browse Source
master
Per A Brodtkorb 8 years ago
parent d16d175190
commit 36f9b69a69
4 changed files with 52 additions and 42 deletions
Show Stats Download Patch File Download Diff File

3
wafo/misc.py
Unescape Escape View File

@ -989,7 +989,8 @@ def findextrema(x):
findcross findcross
crossdef crossdef
''' '''
return findcross(diff(x), 0.0) + 1 x1 = np.atleast_1d(x)
return findcross(diff(x1), 0.0) + 1
def findpeaks(data, n=2, min_h=None, min_p=0.0): def findpeaks(data, n=2, min_h=None, min_p=0.0):

2
wafo/objects.py
Unescape Escape View File

@ -2491,7 +2491,7 @@ class TimeSeries(PlotData):
plt.title('Surface elevation from mean water level (MWL).') plt.title('Surface elevation from mean water level (MWL).')
for ix in range(nsub): for ix in range(nsub):
if nsub > 1: if nsub > 1:
subplot(nsub, 1, ix) subplot(nsub, 1, ix+1)
h_scale = array([tn[ind[0]], tn[ind[-1]]]) h_scale = array([tn[ind[0]], tn[ind[-1]]])
ind2 = where((h_scale[0] <= tn2) & (tn2 <= h_scale[1]))[0] ind2 = where((h_scale[0] <= tn2) & (tn2 <= h_scale[1]))[0]
plot(tn[ind] * dT, xn[ind], sym1) plot(tn[ind] * dT, xn[ind], sym1)

63
wafo/spectrum/core.py
Unescape Escape View File

@ -14,30 +14,30 @@ from scipy.special import erf
from scipy.linalg import toeplitz from scipy.linalg import toeplitz
import scipy.interpolate as interpolate import scipy.interpolate as interpolate
from scipy.interpolate.interpolate import interp1d, interp2d from scipy.interpolate.interpolate import interp1d, interp2d
from ..objects import TimeSeries, mat2timeseries from wafo.objects import TimeSeries, mat2timeseries
from ..interpolate import stineman_interp from wafo.interpolate import stineman_interp
from ..wave_theory.dispersion_relation import w2k # , k2w from wafo.wave_theory.dispersion_relation import w2k # , k2w
from ..containers import PlotData, now from wafo.containers import PlotData, now
from ..misc import sub_dict_select, nextpow2, discretize, JITImport from wafo.misc import (sub_dict_select, nextpow2, discretize, JITImport,
from ..misc import meshgrid, gravity, cart2polar, polar2cart meshgrid, cart2polar, polar2cart, gravity as _gravity)
from ..markov import mctp2rfc, mctp2tc from wafo.markov import mctp2rfc, mctp2tc
from ..kdetools import qlevels from wafo.kdetools import qlevels
# from wafo.transform import TrData # from wafo.transform import TrData
from ..transform.models import TrLinear from wafo.transform.models import TrLinear
from ..plotbackend import plotbackend from wafo.plotbackend import plotbackend
try: try:
from ..gaussian import Rind from wafo.gaussian import Rind
except ImportError: except ImportError:
Rind = None Rind = None
try: try:
from .. import c_library from wafo import c_library
except ImportError: except ImportError:
warnings.warn('Compile the c_library.pyd again!') warnings.warn('Compile the c_library.pyd again!')
c_library = None c_library = None
try: try:
from .. import cov2mod from wafo import cov2mod
except ImportError: except ImportError:
warnings.warn('Compile the cov2mod.pyd again!') warnings.warn('Compile the cov2mod.pyd again!')
cov2mod = None cov2mod = None
@ -102,7 +102,7 @@ def qtf(w, h=inf, g=9.81):
k1p2 = (k_1 + k_2) k1p2 = (k_1 + k_2)
k1m2 = abs(k_1 - k_2) k1m2 = abs(k_1 - k_2)
if 0: # Langley if False: # Langley
p_1 = (-2 * w1p2 * (k12 * g ** 2. - w12 ** 2.) + p_1 = (-2 * w1p2 * (k12 * g ** 2. - w12 ** 2.) +
w_1 * (w_2 ** 4. - g ** 2 * k_2 ** 2) + w_1 * (w_2 ** 4. - g ** 2 * k_2 ** 2) +
w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12) w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12)
@ -914,8 +914,8 @@ class SpecData1D(PlotData):
# ftype = 'w' # ftype = 'w'
# S = ttspec(S,ftype) # S = ttspec(S,ftype)
# end # end
Hm0 = self.characteristic('Hm0') Hm0 = self.characteristic('Hm0')[0]
Tm02 = self.characteristic('Tm02') Tm02 = self.characteristic('Tm02')[0]
if iseed is not None: if iseed is not None:
_set_seed(iseed) # set the the seed _set_seed(iseed) # set the the seed
@ -937,8 +937,7 @@ class SpecData1D(PlotData):
num_waves = 10000 # Typical number of waves in 30 hour seastate num_waves = 10000 # Typical number of waves in 30 hour seastate
Amax = sqrt(2 * log(num_waves)) * Hm0 / 4 Amax = sqrt(2 * log(num_waves)) * Hm0 / 4
fLimitLo = sqrt( fLimitLo = sqrt(gravity * tanh(kbar * water_depth) * Amax / water_depth ** 3)
gravity * tanh(kbar * water_depth) * Amax / water_depth ** 3)
freq = S.args freq = S.args
eps = finfo(float).eps eps = finfo(float).eps
@ -955,7 +954,7 @@ class SpecData1D(PlotData):
# Fs = 2*freq(end)+eps # sampling frequency # Fs = 2*freq(end)+eps # sampling frequency
for ix in range(max_sim): for ix in range(max_sim):
x2, x1 = self.sim_nl(ns=np, cases=cases, dt=None, iseed=iseed, x2, x1 = self.sim_nl(ns=ns, cases=cases, dt=None, iseed=iseed,
method=method, fnlimit=fn_limit, method=method, fnlimit=fn_limit,
output='timeseries') output='timeseries')
x2.data -= x1.data # x2(:,2:end) = x2(:,2:end) -x1(:,2:end) x2.data -= x1.data # x2(:,2:end) = x2(:,2:end) -x1(:,2:end)
@ -968,10 +967,12 @@ class SpecData1D(PlotData):
# %[tf21,fi] = tfe(x2(:,2),x1(:,2),1024,Fs,[],512) # %[tf21,fi] = tfe(x2(:,2),x1(:,2),1024,Fs,[],512)
# %Hw11 = interp1q(fi,tf21.*conj(tf21),freq) # %Hw11 = interp1q(fi,tf21.*conj(tf21),freq)
if True: if True:
Hw1 = exp(interp1d(log(abs(S1.data / S2.data)), S2.args)(freq)) Hw1 = exp(interp1d(log(abs(S1.data / (S2.data + 1e-5))), S2.args,
fill_value=0, bounds_error=False)(freq))
else: else:
# Geometric mean # Geometric mean
fun = interp1d(log(abs(S1.data / S2.data)), S2.args) fun = interp1d(log(abs(S1.data / (S2.data + 1e-5))), S2.args,
fill_value=0, bounds_error=False)
Hw1 = exp((fun(freq) + log(Hw2)) / 2) Hw1 = exp((fun(freq) + log(Hw2)) / 2)
# end # end
# Hw1 = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2 # Hw1 = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2
@ -1011,6 +1012,8 @@ class SpecData1D(PlotData):
plotbackend.figure(2) plotbackend.figure(2)
plotbackend.semilogy(freq, abs(Hw12), 'b') plotbackend.semilogy(freq, abs(Hw12), 'b')
plotbackend.title('Hw-HwOld') plotbackend.title('Hw-HwOld')
plotbackend.show('hold')
# figtile # figtile
# end # end
@ -3068,9 +3071,9 @@ class SpecData1D(PlotData):
>>> S = S0.copy() >>> S = S0.copy()
>>> me, va, sk, ku = S.stats_nl(moments='mvsk') >>> me, va, sk, ku = S.stats_nl(moments='mvsk')
>>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, >>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4)
... ysigma=Hs/4)
>>> ys = wo.mat2timeseries(S.sim(ns=2**13)) >>> ys = wo.mat2timeseries(S.sim(ns=2**13))
>>> ys
>>> g0, gemp = ys.trdata() >>> g0, gemp = ys.trdata()
>>> t0 = g0.dist2gauss() >>> t0 = g0.dist2gauss()
>>> t1 = S0.testgaussian(ns=2**13, cases=50) >>> t1 = S0.testgaussian(ns=2**13, cases=50)
@ -3363,6 +3366,11 @@ class SpecData1D(PlotData):
self.data = intfun(self.args) self.data = intfun(self.args)
self.data = self.data.clip(0) # clip negative values to 0 self.data = self.data.clip(0) # clip negative values to 0
def interp(self, dt):
S = self.copy()
S.resample(dt)
return S
def normalize(self, gravity=9.81): def normalize(self, gravity=9.81):
''' '''
Normalize a spectral density such that m0=m2=1 Normalize a spectral density such that m0=m2=1
@ -4036,7 +4044,7 @@ class SpecData2D(PlotData):
if nr > 0: if nr > 0:
vec = [] vec = []
g = np.atleast_1d(S1.__dict__.get('g', gravity())) g = np.atleast_1d(S1.__dict__.get('g', _gravity()))
# maybe different normalization in x and y => diff. g # maybe different normalization in x and y => diff. g
kx = w ** 2 / g[0] kx = w ** 2 / g[0]
ky = w ** 2 / g[-1] ky = w ** 2 / g[-1]
@ -4231,7 +4239,8 @@ def test_mm_pdf():
S = sm.SpecData1D(Sj(w), w) # Alternatively do it manually S = sm.SpecData1D(Sj(w), w) # Alternatively do it manually
S0 = S.to_linspec() S0 = S.to_linspec()
mm = S.to_mm_pdf() mm = S.to_mm_pdf()
mm.plot()
plotbackend.show()
def test_docstrings(): def test_docstrings():
import doctest import doctest
@ -4239,6 +4248,6 @@ def test_docstrings():
if __name__ == '__main__': if __name__ == '__main__':
test_docstrings() # test_docstrings()
# test_mm_pdf() test_mm_pdf()
# main() # main()

24
wafo/spectrum/models.py
Unescape Escape View File

@ -36,9 +36,9 @@ import scipy.integrate as integrate
import scipy.special as sp import scipy.special as sp
from scipy.fftpack import fft from scipy.fftpack import fft
import numpy as np import numpy as np
from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, from numpy import (inf, atleast_1d, newaxis, minimum, maximum, array,
asarray, exp, log, sqrt, where, pi, arange, linspace, sin, asarray, exp, log, sqrt, where, pi, arange, linspace, sin,
cos, abs, sinh, isfinite, mod, expm1, tanh, cosh, finfo, cos, sinh, isfinite, mod, expm1, tanh, cosh, finfo,
ones, ones_like, isnan, zeros_like, flatnonzero, sinc, ones, ones_like, isnan, zeros_like, flatnonzero, sinc,
hstack, vstack, real, flipud, clip) hstack, vstack, real, flipud, clip)
from ..wave_theory.dispersion_relation import w2k, k2w # @UnusedImport from ..wave_theory.dispersion_relation import w2k, k2w # @UnusedImport
@ -498,7 +498,7 @@ class Jonswap(ModelSpectrum):
>>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1) >>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
>>> S2 = wsm.Bretschneider(Hm0=7, Tp=11) >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
>>> w = plb.linspace(0,5) >>> w = plb.linspace(0,5)
>>> all(abs(S(w)-S2(w))<1.e-7) >>> all(np.abs(S(w)-S2(w))<1.e-7)
True True
h = plb.plot(w,S(w)) h = plb.plot(w,S(w))
@ -1313,7 +1313,7 @@ class Wallop(Bretschneider):
wp = 2. * pi / Tp wp = 2. * pi / Tp
kp = w2k(wp, 0, inf)[0] # wavenumber at peak frequency kp = w2k(wp, 0, inf)[0] # wavenumber at peak frequency
Lp = 2. * pi / kp # wave length at the peak frequency Lp = 2. * pi / kp # wave length at the peak frequency
N = abs((log(2. * pi ** 2.) + 2 * log(Hm0 / 4) - N = np.abs((log(2. * pi ** 2.) + 2 * log(Hm0 / 4) -
2.0 * log(Lp)) / log(2)) 2.0 * log(Lp)) / log(2))
super(Wallop, self).__init__(Hm0, Tp, N, M, chk_seastate) super(Wallop, self).__init__(Hm0, Tp, N, M, chk_seastate)
@ -1796,7 +1796,7 @@ class Spreading(object):
""" Returns the solution of r1 = x. """ Returns the solution of r1 = x.
""" """
X = r1 X = r1
if any(X >= 1): if np.any(X >= 1):
raise ValueError('POISSON spreading: X value must be less than 1') raise ValueError('POISSON spreading: X value must be less than 1')
return X return X
@ -1822,9 +1822,9 @@ class Spreading(object):
A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0) A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0)
ix = flatnonzero( ix = flatnonzero(
(abs(da) > sqrt(_EPS) * abs(A)) * (abs(da) > sqrt(_EPS))) (np.abs(da) > sqrt(_EPS) * np.abs(A)) * (np.abs(da) > sqrt(_EPS)))
if ix.size == 0: if ix.size == 0:
if any(A > pi): if np.any(A > pi):
raise ValueError( raise ValueError(
'BOX-CAR spreading: The A value must be less than pi') 'BOX-CAR spreading: The A value must be less than pi')
return A.clip(min=1e-16, max=pi) return A.clip(min=1e-16, max=pi)
@ -1872,7 +1872,7 @@ class Spreading(object):
def fun(x): def fun(x):
return 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix] return 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix]
for ix in ix0: for ix in ix0:
B[ix] = abs(optimize.fsolve(fun, B0[ix])) B[ix] = np.abs(optimize.fsolve(fun, B0[ix]))
return B return B
def fourier2d(self, r1): def fourier2d(self, r1):
@ -1953,14 +1953,14 @@ class Spreading(object):
self._frequency_independent_spread) self._frequency_independent_spread)
s = spread(wn) s = spread(wn)
if any(s < 0): if np.any(s < 0):
raise ValueError('The COS2S spread parameter, S(w), ' + raise ValueError('The COS2S spread parameter, S(w), ' +
'value must be larger than 0') 'value must be larger than 0')
if self.type[0] == 'c': # cos2s if self.type[0] == 'c': # cos2s
s_par = s s_par = s
else: else:
# First Fourier coefficient of the directional spreading function. # First Fourier coefficient of the directional spreading function.
r1 = abs(s / (s + 1)) r1 = np.abs(s / (s + 1))
# Find distribution parameter from first Fourier coefficient. # Find distribution parameter from first Fourier coefficient.
s_par = self.fourier2distpar(r1) s_par = self.fourier2distpar(r1)
if self.method is not None: if self.method is not None:
@ -1986,8 +1986,8 @@ class Spreading(object):
@staticmethod @staticmethod
def _check_theta(theta): def _check_theta(theta):
L = abs(theta[-1] - theta[0]) L = np.abs(theta[-1] - theta[0])
if abs(L - 2 * np.pi) > _EPS: if np.abs(L - 2 * np.pi) > _EPS:
raise ValueError('theta must cover all angles -pi -> pi') raise ValueError('theta must cover all angles -pi -> pi')
nt = len(theta) nt = len(theta)
if nt < 40: if nt < 40:

Write Preview
Loading…
Cancel
Save