From 9522c6c24fcb7c23974a07c19b35e7202a41aaa6 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Sun, 14 Feb 2016 15:10:49 +0100 Subject: [PATCH] Cleanup of code --- wafo/covariance/core.py | 8 +- wafo/covariance/tests/__init__.py | 1 - wafo/gaussian.py | 12 +- wafo/integrate.py | 232 +----- wafo/spectrum/__init__.py | 2 +- wafo/spectrum/core.py | 1008 ++++++++++-------------- wafo/spectrum/models.py | 91 +-- wafo/spectrum/tests/__init__.py | 1 - wafo/spectrum/tests/test_specdata1d.py | 372 ++++----- wafo/stats/_distn_infrastructure.py | 23 +- wafo/tests/__init__.py | 1 - wafo/tests/conftest.py | 2 +- wafo/tests/test_integrate.py | 63 +- wafo/transform/models.py | 4 +- wafo/wave_theory/core.py | 174 ++-- 15 files changed, 831 insertions(+), 1163 deletions(-) diff --git a/wafo/covariance/core.py b/wafo/covariance/core.py index abd6fe3..afb88b2 100644 --- a/wafo/covariance/core.py +++ b/wafo/covariance/core.py @@ -14,7 +14,7 @@ note : Memorandum string. date : Date and time of creation or change. ''' -from __future__ import division +from __future__ import division, absolute_import import warnings import numpy as np from numpy import (zeros, ones, sqrt, inf, where, nan, @@ -27,9 +27,9 @@ from scipy.linalg import toeplitz, lstsq from scipy import sparse from pylab import stineman_interp -from wafo.containers import PlotData -from wafo.misc import sub_dict_select, nextpow2 # , JITImport -import wafo.spectrum as _wafospec +from ..containers import PlotData +from ..misc import sub_dict_select, nextpow2 # , JITImport +from .. import spectrum as _wafospec from scipy.sparse.linalg.dsolve.linsolve import spsolve from scipy.sparse.base import issparse from scipy.signal.windows import parzen diff --git a/wafo/covariance/tests/__init__.py b/wafo/covariance/tests/__init__.py index 8b13789..e69de29 100644 --- a/wafo/covariance/tests/__init__.py +++ b/wafo/covariance/tests/__init__.py @@ -1 +0,0 @@ - diff --git a/wafo/gaussian.py b/wafo/gaussian.py index ae4dbd1..2a3a896 100644 --- a/wafo/gaussian.py +++ b/wafo/gaussian.py @@ -355,8 +355,8 @@ class Rind(object): dev = sqrt(diag(BIG)) # std ind = nonzero(indI[1:] > -1)[0] infin = repeat(2, len(indI) - 1) - infin[ind] = (2 - (Bup[0, ind] > infinity * dev[indI[ind + 1]]) - - 2 * (Blo[0, ind] < -infinity * dev[indI[ind + 1]])) + infin[ind] = (2 - (Bup[0, ind] > infinity * dev[indI[ind + 1]]) - + 2 * (Blo[0, ind] < -infinity * dev[indI[ind + 1]])) Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]]) Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]]) @@ -992,10 +992,10 @@ def prbnorm2d(a, b, r): infin = np.repeat(2, 2) - (upper > infinity) - 2 * (lower < -infinity) if np.all(infin == 2): - return (bvd(lower[0], lower[1], correl) - - bvd(upper[0], lower[1], correl) - - bvd(lower[0], upper[1], correl) - + bvd(upper[0], upper[1], correl)) + return (bvd(lower[0], lower[1], correl) - + bvd(upper[0], lower[1], correl) - + bvd(lower[0], upper[1], correl) + + bvd(upper[0], upper[1], correl)) elif (infin[0] == 2 and infin[1] == 1): return (bvd(lower[0], lower[1], correl) - bvd(upper[0], lower[1], correl)) diff --git a/wafo/integrate.py b/wafo/integrate.py index 5783bab..99ccc8d 100644 --- a/wafo/integrate.py +++ b/wafo/integrate.py @@ -266,8 +266,8 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3): fp[i] = 4 * fp[i - 1] # Richardson extrapolation for k in range(i): - rom[two, k + 1] = rom[two, k] + \ - (rom[two, k] - rom[one, k]) / (fp[k] - 1) + rom[two, k + 1] = (rom[two, k] + + (rom[two, k] - rom[one, k]) / (fp[k] - 1)) Ih1 = Ih2 Ih2 = Ih4 @@ -1119,6 +1119,8 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17): ''' # Author: jonas.lundgren@saabgroup.com, 2009. license BSD # Order limits (required if infinite limits) + a = np.asarray(a) + b = np.asarray(b) if a == b: Q = b - a err = b - a @@ -1138,17 +1140,17 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17): # Change of variable if np.isfinite(a) & np.isinf(b): # a to inf - fun1 = lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2 - [Q, err] = quadgr(fun1, 0, 1, abseps) + [Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2, + 0, 1, abseps) elif np.isinf(a) & np.isfinite(b): # -inf to b - fun2 = lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2 - [Q, err] = quadgr(fun2, -1, 0, abseps) + [Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2, + -1, 0, abseps) else: # -inf to inf - fun1 = lambda t: fun(t / (1 - t)) / (1 - t) ** 2 - fun2 = lambda t: fun(t / (1 + t)) / (1 + t) ** 2 - [Q1, err1] = quadgr(fun1, 0, 1, abseps / 2) - [Q2, err2] = quadgr(fun2, -1, 0, abseps / 2) + [Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2, + 0, 1, abseps / 2) + [Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2, + -1, 0, abseps / 2) Q = Q1 + Q2 err = err1 + err2 @@ -1170,9 +1172,9 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17): dtype = np.result_type(fun(a), fun(b)) # Initiate vectors - Q0 = zeros(max_iter, dtype=dtype) # Quadrature - Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation - Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation + Q0 = zeros(max_iter, dtype=dtype) # Quadrature + Q1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation + Q2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation # One interval hh = (b - a) / 2 # Half interval length @@ -1187,8 +1189,8 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17): hh = hh / 2 x = np.hstack([x + a, x + b]) / 2 # Quadrature - Q0[k] = hh * \ - np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0), axis=0) + Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0), + axis=0) # Richardson extrapolation if k >= 5: @@ -1424,207 +1426,7 @@ def test_docstrings(): doctest.testmod() -# def levin_integrate(): -# ''' An oscillatory integral -# Sheehan Olver, December 2010 -# -# -# (Chebfun example quad/LevinIntegrate.m) -# -# This example computes the highly oscillatory integral of -# -# f * exp( 1i * w * g ), -# -# over (0,1) using the Levin method [1]. This method computes the integral -# by rewriting it as an ODE -# -# u' + 1i * w * g' u = f, -# -# so that the indefinite integral of f * exp( 1i * w * g ) is -# -# u * exp( 1i * w * g ). -# -# -# -# We use as an example -# -# f = 1 / ( x + 2 ); -# g = cos( x - 2 ); -# w = 100000; -# -# # -# References: -# -# [1] Levin, D., Procedures for computing one and two-dimensional integrals -# of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538 -# ''' -# exp = np.exp -# domain=[0, 1] -# x = Chebfun.identity(domain=domain) -# f = 1./(x+2) -# g = np.cos(x-2) -# D = np.diff(domain) -# -# -# # Here is are plots of this integrand, with w = 100, in complex space -# w = 100; -# line_opts = dict(line_width=1.6) -# font_opts = dict(font_size= 14) -# # -# -# intg = f*exp(1j*w*g) -# xs, ys, xi, yi, d = intg.plot_data(1000) -# #intg.plot(with_interpolation_points=True) -# #xi = np.linspace(0, 1, 1024) -# # plt.plot(xs, ys) # , **line_opts) -# # plt.plot(xi, yi, 'r.') -# # #axis equal -# # plt.title('Complex plot of integrand') #,**font_opts) -# # plt.show('hold') -# ## -# # and of just the real part -# # intgr = np.real(intg) -# # xs, ys, xi, yi, d = intgr.plot_data(1000) -# #intgr.plot() -# # plt.plot(xs, np.real(ys)) # , **line_opts) -# # plt.plot(xi, np.real(yi), 'r.') -# #axis equal -# # plt.title('Real part of integrand') #,**font_opts) -# # plt.show('hold') -# -# ## -# # The Levin method will be accurate for large and small w, and the time -# # taken is independent of w. Here we take a reasonably large value of w. -# w = 1000; -# intg = f*exp(1j*w*g) -# val0 = np.sum(intg) -# # val1 = sum(intg) -# print(val0) -# ## -# # Start timing -# #tic -# -# ## -# # Construct the operator L -# L = D + 1j*w*np.diag(g.differentiate()) -# -# ## -# # From asymptotic analysis, we know that there exists a solution to the -# # equation which is non-oscillatory, though we do not know what initial -# # condition it satisfies. Thus we find a particular solution to this -# # equation with no boundary conditions. -# -# u = L / f -# -# ## -# # Because L is a differential operator with derivative order 1, \ expects -# # it to be given a boundary condition, which is why the warning message is -# # displayed. However, this doesn't cause any problems: though there are, -# # in fact, a family of solutions to the ODE without boundary conditions -# # due to the kernel -# # -# # exp(- 1i * w * g), -# # -# # it does not actually matter which particular solution is computed. -# # Non-uniqueness is also not an issue: \ in matlab is least squares, hence -# # does not require uniqueness. The existence of a non-oscillatory solution -# # ensures that \ converges to a u with length independent of w. -# # -# # One could prevent the warning by applying a boundary condition consistent -# # with the rest of the system, that is -# # L.lbc = {L(1,:),f(0)}; -# -# ## -# # Now we evaluate the antiderivative at the endpoints to obtain the -# # integral. -# -# u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0)) -# -# #toc -# -# -# ## -# # Here is a way to compute the integral using Clenshaw--Curtis quadrature. -# # As w becomes large, this takes an increasingly long time as the -# # oscillations must be resolved. -# -# #tic -# sum( f*exp(1j*w*g) ) -# #toc - - -# aLevinTQ[omega_,a_,b_,f_,g_,nu_,wprec_,prm_,test_,basis_,np_]:= -# -# Module[{r,betam,A,AA,BB,S,F,w=N[omega, wprec]}, -# M=Length[nu]-1; -# PB[k_,t_]:=If[basis==1,t^k,ChebyshevT[k,t]]; -# -# ff[t_]:=((b-a)/2)*f[(b-a)*t/2+(a+b)/2]; -# -# gg[t_]:=g[(b-a)*t/2+(a+b)/2]; -# dgg[t_]:=Derivative[1][gg][t]; -# -# If[test==0, betam=Min[Abs[dgg[-1]*w], Abs[dgg[1]*w]]; -# While[prm*M/betam >=1, betam=2*betam]]; -# If[test>0,x[k_]:=N[Cos[k*Pi/M], wprec],x[k_]:= -# Which[kprm*M, N[1-(M-k)/betam, wprec]]]; -# -# Psi[k_,t_]:=Derivative[0,1][PB][k,t]+I*w*dgg[t]*PB[k,t]; -# -# ne[j_]:=nu[[j+1]]; S[-1]=0; S[j_]:=Sum[ne[i],{i,0,j}]; -# nn=S[M]-1; -# A=ConstantArray[0,{nn+1,nn+1}]; -# F=ConstantArray[0,nn+1]; r=0; -# While[rx[r]], -# {k,0,S[M]-1}],{j,S[r-1],S[r]-1}]; -# -# Do[BB[j]=Limit[Derivative[Mod[j-S[r-1],ne[r]]][ff][t], -# t->x[r]],{j,S[r-1],S[r]-1}]; -# Do[F[[j]]=BB[j-1],{j,S[r-1]+1,S[r]}]; -# Do[Do[A[[j,k]]=AA[j-1,k-1],{k,1,S[M]}],{j,S[r-1]+1,S[r]}]; -# r=r+1;]; (*sv=SingularValueList[N[A,wprec]]; -# con=sv[[1]]/sv[[-1]]; Print["cond2(A)= ",N[con,3]];*) -# LS=Block[{MaxExtraPrecision=0},LeastSquares[N[A, wprec],F]]; -# vvv[t_]:=Sum[LS[[k+1]]*PB[k,t], {k,0,nn}]; -# NR=vvv[1]*Exp[I*w*gg[1]]-vvv[-1]*Exp[I*w*gg[-1]]; -# Print["n=", np+ii+2s-2, ", Result= ", N[NR, wprec/2+5]]; -# If[ii==0,PR=NR];]; -# (* End of subroutine aLevinTQ /A.I. Hascelik, July 2013/ *) -# -# def main_levin(): -# a=1; b=2; -# omega=100; -# prm=1/2; -# f[t_]=Exp[4t] -# g[t_]=t+Exp[4t]*Gamma[t] -# -# dg[t_]:=Derivative[1][g][t]; -# -# prec=16 -# wprec=2*prec -# delta = min(abs(omega*dg(a)), abs(omega*dg(b))) -# alpha = min(abs(omega*g(a)), abs(omega*g(b))) -# s=1; #(*if s>1, the integral is computed by Q_s^L*) -# test= 1 if delta<10 or alpha <=10 or s>1 else 0 -# -# m = 1 if s>1 else np.floor(prec/max(np.log10(beta+1),1)+2) -# nc = 2*m+1 #(*or np=2m, number of collocation points*) -# basis=1; # (*take basis=0 for the Chebysev basis*) -# for ii in range(0, 2, 2): -# nu = np.ones((nc+ii,)) # ConstantArray[1,nc+ii]; -# nu[0] = s -# nu[-1] = s -# #nu[[1]]=s; -# #nu[[-1]]=s; -# aLevinTQ[omega,a,b,f,g,nu,wprec,prm,test,basis,nc], -# #{ii,0,2,2}]; -# Print["Error= ",Abs[NR-PR]]; - - if __name__ == '__main__': - # levin_integrate() test_docstrings() # qdemo(np.exp, 0, 3, plot_error=True) # plt.show('hold') diff --git a/wafo/spectrum/__init__.py b/wafo/spectrum/__init__.py index 8106d0d..a112250 100644 --- a/wafo/spectrum/__init__.py +++ b/wafo/spectrum/__init__.py @@ -7,4 +7,4 @@ Spectrum package in WAFO Toolbox. from __future__ import absolute_import from .core import * from . import models -from wafo.wave_theory import dispersion_relation +from ..wave_theory import dispersion_relation diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index cfdfd91..886f365 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -1,6 +1,4 @@ from __future__ import absolute_import, division -from wafo.misc import meshgrid, gravity, cart2polar, polar2cart -from wafo.objects import TimeSeries, mat2timeseries import warnings import os import numpy as np @@ -15,15 +13,18 @@ from scipy.integrate import simps, trapz from scipy.special import erf from scipy.linalg import toeplitz import scipy.interpolate as interpolate +from scipy.interpolate.interpolate import interp1d, interp2d +from ..misc import meshgrid, gravity, cart2polar, polar2cart +from ..objects import TimeSeries, mat2timeseries from ..interpolate import stineman_interp - from ..wave_theory.dispersion_relation import w2k # , k2w from ..containers import PlotData, now -# , tranproc from ..misc import sub_dict_select, nextpow2, discretize, JITImport -# from wafo.graphutil import cltext from ..kdetools import qlevels -from scipy.interpolate.interpolate import interp1d + +# from wafo.transform import TrData +from ..transform.models import TrLinear +from ..plotbackend import plotbackend try: from ..gaussian import Rind @@ -40,20 +41,16 @@ except ImportError: warnings.warn('Compile the cov2mod.pyd again!') cov2mod = None - -# from wafo.transform import TrData -from ..transform.models import TrLinear -from ..plotbackend import plotbackend - - # Trick to avoid error due to circular import - _WAFOCOV = JITImport('wafo.covariance') __all__ = ['SpecData1D', 'SpecData2D', 'plotspec'] +_EPS = np.finfo(float).eps + + def _set_seed(iseed): '''Set seed of random generator''' if iseed is not None: @@ -157,477 +154,369 @@ def qtf(w, h=inf, g=9.81): def plotspec(specdata, linetype='b-', flag=1): + ''' + PLOTSPEC Plot a spectral density + + Parameters + ---------- + S : SpecData1D or SpecData2D object + defining spectral density. + linetype : string + defining color and linetype, see plot for possibilities + flag : scalar integer + defining the type of plot + 1D: + 1 plots the density, S, (default) + 2 plot 10log10(S) + 3 plots both the above plots + 2D: + Directional spectra: S(w,theta), S(f,theta) + 1 polar plot S (default) + 2 plots spectral density and the directional + spreading, int S(w,theta) dw or int S(f,theta) df + 3 plots spectral density and the directional + spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df + 4 mesh of S + 5 mesh of S in polar coordinates + 6 contour plot of S + 7 filled contour plot of S + Wavenumber spectra: S(k1,k2) + 1 contour plot of S (default) + 2 filled contour plot of S + + Example + ------- + >>> import numpy as np + >>> import wafo.spectrum as ws + >>> Sj = ws.models.Jonswap(Hm0=3, Tp=7) + >>> S = Sj.tospecdata() + >>> ws.plotspec(S,flag=1) + + S = demospec('dir'); S2 = mkdspec(jonswap,spreading); + plotspec(S,2), hold on + # Same as previous fig. due to frequency independent spreading + plotspec(S,3,'g') + # Not the same as previous figs. due to frequency dependent spreading + plotspec(S2,2,'r') + plotspec(S2,3,'m') + # transform from angular frequency and radians to frequency and degrees + Sf = ttspec(S,'f','d'); clf + plotspec(Sf,2), + + See also + dat2spec, createspec, simpson + ''' pass -# ''' -# PLOTSPEC Plot a spectral density -# -# Parameters -# ---------- -# S : SpecData1D or SpecData2D object -# defining spectral density. -# linetype : string -# defining color and linetype, see plot for possibilities -# flag : scalar integer -# defining the type of plot -# 1D: -# 1 plots the density, S, (default) -# 2 plot 10log10(S) -# 3 plots both the above plots -# 2D: -# Directional spectra: S(w,theta), S(f,theta) -# 1 polar plot S (default) -# 2 plots spectral density and the directional -# spreading, int S(w,theta) dw or int S(f,theta) df -# 3 plots spectral density and the directional -# spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df -# 4 mesh of S -# 5 mesh of S in polar coordinates -# 6 contour plot of S -# 7 filled contour plot of S -# Wavenumber spectra: S(k1,k2) -# 1 contour plot of S (default) -# 2 filled contour plot of S +# # label the contour levels +# txtFlag = 0 +# LegendOn = 1 # -# Example -# ------- -# >>> import numpy as np -# >>> import wafo.spectrum as ws -# >>> Sj = ws.models.Jonswap(Hm0=3, Tp=7) -# >>> S = Sj.tospecdata() -# >>> ws.plotspec(S,flag=1) +# ftype = specdata.freqtype # options are 'f' and 'w' and 'k' +# data = specdata.data +# if data.ndim == 2: +# freq = specdata.args[1] +# theta = specdata.args[0] +# else: +# freq = specdata.args +# if isinstance(specdata.args, (list, tuple)): # -# S = demospec('dir'); S2 = mkdspec(jonswap,spreading); -# plotspec(S,2), hold on -# # Same as previous fig. due to frequency independent spreading -# plotspec(S,3,'g') -# # Not the same as previous figs. due to frequency dependent spreading -# plotspec(S2,2,'r') -# plotspec(S2,3,'m') -# % transform from angular frequency and radians to frequency and degrees -# Sf = ttspec(S,'f','d'); clf -# plotspec(Sf,2), +# if ftype == 'w': +# xlbl_txt = 'Frequency [rad/s]' +# ylbl1_txt = 'S(w) [m^2 s / rad]' +# ylbl3_txt = 'Directional Spectrum' +# zlbl_txt = 'S(w,\theta) [m^2 s / rad^2]' +# funit = ' [rad/s]' +# Sunit = ' [m^2 s / rad]' +# elif ftype == 'f': +# xlbl_txt = 'Frequency [Hz]' +# ylbl1_txt = 'S(f) [m^2 s]' +# ylbl3_txt = 'Directional Spectrum' +# zlbl_txt = 'S(f,\theta) [m^2 s / rad]' +# funit = ' [Hz]' +# Sunit = ' [m^2 s ]' +# elif ftype == 'k': +# xlbl_txt = 'Wave number [rad/m]' +# ylbl1_txt = 'S(k) [m^3/ rad]' +# funit = ' [rad/m]' +# Sunit = ' [m^3 / rad]' +# ylbl4_txt = 'Wave Number Spectrum' # -# See also dat2spec, createspec, simpson -# ''' +# else: +# raise ValueError('Frequency type unknown') # -# label the contour levels -# txtFlag = 0 -# LegendOn = 1 # +# if hasattr(specdata, 'norm') and specdata.norm : +# Sunit=[] +# funit = [] +# ylbl1_txt = 'Normalized Spectral density' +# ylbl3_txt = 'Normalized Directional Spectrum' +# ylbl4_txt = 'Normalized Wave Number Spectrum' +# if ftype == 'k': +# xlbl_txt = 'Normalized Wave number' +# else: +# xlbl_txt = 'Normalized Frequency' # -# ftype = specdata.freqtype #options are 'f' and 'w' and 'k' -# data = specdata.data -# if data.ndim == 2: -# freq = specdata.args[1] -# theta = specdata.args[0] -# else: -# freq = specdata.args -# if isinstance(specdata.args, (list, tuple)): +# ylbl2_txt = 'Power spectrum (dB)' # -# if ftype == 'w': -# xlbl_txt = 'Frequency [rad/s]' -# ylbl1_txt = 'S(w) [m^2 s / rad]' -# ylbl3_txt = 'Directional Spectrum' -# zlbl_txt = 'S(w,\theta) [m^2 s / rad^2]' -# funit = ' [rad/s]' -# Sunit = ' [m^2 s / rad]' -# elif ftype == 'f': -# xlbl_txt = 'Frequency [Hz]' -# ylbl1_txt = 'S(f) [m^2 s]' -# ylbl3_txt = 'Directional Spectrum' -# zlbl_txt = 'S(f,\theta) [m^2 s / rad]' -# funit = ' [Hz]' -# Sunit = ' [m^2 s ]' -# elif ftype == 'k': -# xlbl_txt = 'Wave number [rad/m]' -# ylbl1_txt = 'S(k) [m^3/ rad]' -# funit = ' [rad/m]' -# Sunit = ' [m^3 / rad]' -# ylbl4_txt = 'Wave Number Spectrum' +# phi = specdata.phi # -# else: -# raise ValueError('Frequency type unknown') +# spectype = specdata.type.lower() +# stype = spectype[-3::] +# if stype in ('enc', 'req', 'k1d') : #1D plot +# Fn = freq[-1] # Nyquist frequency +# indm = findpeaks(data, n=4) +# maxS = data.max() +# if isfield(S,'CI') && ~isempty(S.CI): +# maxS = maxS*S.CI(2) +# txtCI = [num2str(100*S.p), '% CI'] +# #end # +# Fp = freq[indm]# %peak frequency/wave number # -# if hasattr(specdata, 'norm') and specdata.norm : -# Sunit=[] -# funit = [] -# ylbl1_txt = 'Normalized Spectral density' -# ylbl3_txt = 'Normalized Directional Spectrum' -# ylbl4_txt = 'Normalized Wave Number Spectrum' -# if ftype == 'k': -# xlbl_txt = 'Normalized Wave number' -# else: -# xlbl_txt = 'Normalized Frequency' +# if len(indm) == 1: +# txt = [('fp = %0.2g' % Fp) + funit] +# else: +# txt = [] +# for i, fp in enumerate(Fp.tolist()): +# txt.append(('fp%d = %0.2g' % (i, fp)) + funit) # -# ylbl2_txt = 'Power spectrum (dB)' +# txt = ''.join(txt) +# if (flag == 3): +# plotbackend.subplot(2, 1, 1) +# if (flag == 1) or (flag == 3):# Plot in normal scale +# plotbackend.plot(np.vstack([Fp, Fp]), +# np.vstack([zeros(len(indm)), data.take(indm)]), +# ':', label=txt) +# plotbackend.plot(freq, data, linetype) +# specdata.labels.labelfig() +# if isfield(S,'CI'): +# plot(freq,S.S*S.CI(1), 'r:' ) +# plot(freq,S.S*S.CI(2), 'r:' ) # -# phi = specdata.phi +# a = plotbackend.axis() # -# spectype = specdata.type.lower() -# stype = spectype[-3::] -# if stype in ('enc', 'req', 'k1d') : #1D plot -# Fn = freq[-1] # Nyquist frequency -# indm = findpeaks(data, n=4) -# maxS = data.max() -# if isfield(S,'CI') && ~isempty(S.CI), -# maxS = maxS*S.CI(2) -# txtCI = [num2str(100*S.p), '% CI'] -# end +# a1 = Fn +# if (Fp > 0): +# a1 = max(min(Fn, 10 * max(Fp)), a[1]) # -# Fp = freq[indm]# %peak frequency/wave number +# plotbackend.axis([0, a1 , 0, max(1.01 * maxS, a[3])]) +# plotbackend.title('Spectral density') +# plotbackend.xlabel(xlbl_txt) +# plotbackend.ylabel(ylbl1_txt) # -# if len(indm) == 1: -# txt = [('fp = %0.2g' % Fp) + funit] -# else: -# txt = [] -# for i, fp in enumerate(Fp.tolist()): -# txt.append(('fp%d = %0.2g' % (i, fp)) + funit) # -# txt = ''.join(txt) -# if (flag == 3): -# plotbackend.subplot(2, 1, 1) -# if (flag == 1) or (flag == 3):# Plot in normal scale -# plotbackend.plot(np.vstack([Fp, Fp]), -# np.vstack([zeros(len(indm)), data.take(indm)]), -# ':', label=txt) -# plotbackend.plot(freq, data, linetype) -# specdata.labels.labelfig() -# if isfield(S,'CI'), -# plot(freq,S.S*S.CI(1), 'r:' ) -# plot(freq,S.S*S.CI(2), 'r:' ) +# if (flag == 3): +# plotbackend.subplot(2, 1, 2) # -# a = plotbackend.axis() +# if (flag == 2) or (flag == 3) : # Plot in logaritmic scale +# ind = np.flatnonzero(data > 0) # -# a1 = Fn -# if (Fp > 0): -# a1 = max(min(Fn, 10 * max(Fp)), a[1]) +# plotbackend.plot(np.vstack([Fp, Fp]), +# np.vstack((min(10 * log10(data.take(ind) / +# maxS)).repeat(len(Fp)), +# 10 * log10(data.take(indm) / maxS))), ':',label=txt) +# hold on +# if isfield(S,'CI'): +# plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' ) +# plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' ) # -# plotbackend.axis([0, a1 , 0, max(1.01 * maxS, a[3])]) -# plotbackend.title('Spectral density') -# plotbackend.xlabel(xlbl_txt) -# plotbackend.ylabel(ylbl1_txt) +# plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype) # +# a = plotbackend.axis() # -# if (flag == 3): -# plotbackend.subplot(2, 1, 2) +# a1 = Fn +# if (Fp > 0): +# a1 = max(min(Fn, 10 * max(Fp)), a[1]) # -# if (flag == 2) or (flag == 3) : # Plot in logaritmic scale -# ind = np.flatnonzero(data > 0) +# plotbackend.axis([0, a1 , -20, max(1.01 * 10 * log10(1), a[3])]) # -# plotbackend.plot(np.vstack([Fp, Fp]), -# np.vstack((min(10 * log10(data.take(ind) / -# maxS)).repeat(len(Fp)), -# 10 * log10(data.take(indm) / maxS))), ':',label=txt) -# hold on -# if isfield(S,'CI'), -# plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' ) -# plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' ) -# end -# plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype) +# specdata.labels.labelfig() +# plotbackend.title('Spectral density') +# plotbackend.xlabel(xlbl_txt) +# plotbackend.ylabel(ylbl2_txt) # -# a = plotbackend.axis() +# if LegendOn: +# plotbackend.legend() +# if isfield(S,'CI'), +# legend(txt{:},txtCI,1) +# else +# legend(txt{:},1) +# end +# end +# case {'k2d'} +# if plotflag==1, +# [c, h] = contour(freq,S.k2,S.S,'b') +# z_level = clevels(c) # -# a1 = Fn -# if (Fp > 0): -# a1 = max(min(Fn, 10 * max(Fp)), a[1]) # -# plotbackend.axis([0, a1 , -20, max(1.01 * 10 * log10(1), a[3])]) +# if txtFlag==1 +# textstart_x=0.05; textstart_y=0.94 +# cltext1(z_level,textstart_x,textstart_y) +# else +# cltext(z_level,0) +# end +# else +# [c,h] = contourf(freq,S.k2,S.S) +# %clabel(c,h), colorbar(c,h) +# fcolorbar(c) % alternative +# end +# rotate(h,[0 0 1],-phi*180/pi) # -# specdata.labels.labelfig() -# plotbackend.title('Spectral density') -# plotbackend.xlabel(xlbl_txt) -# plotbackend.ylabel(ylbl2_txt) # -# if LegendOn: -# plotbackend.legend() -# if isfield(S,'CI'), -# legend(txt{:},txtCI,1) -# else -# legend(txt{:},1) -# end -# end -# case {'k2d'} -# if plotflag==1, -# [c, h] = contour(freq,S.k2,S.S,'b') -# z_level = clevels(c) # +# xlabel(xlbl_txt) +# ylabel(xlbl_txt) +# title(ylbl4_txt) +# # return +# km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)]) +# axis([-km km -km km]) +# hold on +# plot([0 0],[ -km km],':') +# plot([-km km],[0 0],':') +# axis('square') # -# if txtFlag==1 -# textstart_x=0.05; textstart_y=0.94 -# cltext1(z_level,textstart_x,textstart_y) -# else -# cltext(z_level,0) -# end -# else -# [c,h] = contourf(freq,S.k2,S.S) -# %clabel(c,h), colorbar(c,h) -# fcolorbar(c) % alternative -# end -# rotate(h,[0 0 1],-phi*180/pi) # +# # cltext(z_level) +# # axis('square') +# if ~ih, hold off,end +# case {'dir'} +# thmin = S.theta(1)-phi;thmax=S.theta(end)-phi +# if plotflag==1 % polar plot +# if 0, % alternative but then z_level must be chosen beforehand +# h = polar([0 2*pi],[0 freq(end)]) +# delete(h);hold on +# [X,Y]=meshgrid(S.theta,freq) +# [X,Y]=polar2cart(X,Y) +# contour(X,Y,S.S',lintype) +# else +# if (abs(thmax-thmin)<3*pi), % angle given in radians +# theta = S.theta +# else +# theta = S.theta*pi/180 % convert to radians +# phi = phi*pi/180 +# end +# c = contours(theta,freq,S.S')%,Nlevel) % calculate levels +# if isempty(c) +# c = contours(theta,freq,S.S)%,Nlevel); % calculate levels +# end +# [z_level c] = clevels(c); % find contour levels +# h = polar(c(1,:),c(2,:),lintype); +# rotate(h,[0 0 1],-phi*180/pi) +# end +# title(ylbl3_txt) +# % label the contour levels # +# if txtFlag==1 +# textstart_x = -0.1; textstart_y=1.00; +# cltext1(z_level,textstart_x,textstart_y); +# else +# cltext(z_level,0) +# end # -# xlabel(xlbl_txt) -# ylabel(xlbl_txt) -# title(ylbl4_txt) -# %return -# km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)]) -# axis([-km km -km km]) -# hold on -# plot([0 0],[ -km km],':') -# plot([-km km],[0 0],':') -# axis('square') +# elseif (plotflag==2) || (plotflag==3), +# %ih = ishold; # +# subplot(211) # -# %cltext(z_level) -# %axis('square') -# if ~ih, hold off,end -# case {'dir'} -# thmin = S.theta(1)-phi;thmax=S.theta(end)-phi -# if plotflag==1 % polar plot -# if 0, % alternative but then z_level must be chosen beforehand -# h = polar([0 2*pi],[0 freq(end)]) -# delete(h);hold on -# [X,Y]=meshgrid(S.theta,freq) -# [X,Y]=polar2cart(X,Y) -# contour(X,Y,S.S',lintype) -# else -# if (abs(thmax-thmin)<3*pi), % angle given in radians -# theta = S.theta -# else -# theta = S.theta*pi/180 % convert to radians -# phi = phi*pi/180 -# end -# c = contours(theta,freq,S.S')%,Nlevel) % calculate levels -# if isempty(c) -# c = contours(theta,freq,S.S)%,Nlevel); % calculate levels -# end -# [z_level c] = clevels(c); % find contour levels -# h = polar(c(1,:),c(2,:),lintype); -# rotate(h,[0 0 1],-phi*180/pi) -# end -# title(ylbl3_txt) -# % label the contour levels +# if ih, hold on, end # -# if txtFlag==1 -# textstart_x = -0.1; textstart_y=1.00; -# cltext1(z_level,textstart_x,textstart_y); -# else -# cltext(z_level,0) -# end +# Sf = spec2spec(S,'freq'); % frequency spectrum +# plotspec(Sf,1,lintype) # -# elseif (plotflag==2) || (plotflag==3), -# %ih = ishold; +# subplot(212) # -# subplot(211) +# Dtf = S.S; +# [Nt,Nf] = size(S.S); +# Sf = Sf.S(:).'; +# ind = find(Sf); # -# if ih, hold on, end +# if plotflag==3, %Directional distribution D(theta,freq)) +# Dtf(:,ind) = Dtf(:,ind)./Sf(ones(Nt,1),ind); +# end +# Dtheta = simpson(freq,Dtf,2); %Directional spreading, D(theta) +# Dtheta = Dtheta/simpson(S.theta,Dtheta); % make sure int D(theta)dtheta = 1 +# [y,ind] = max(Dtheta); +# Wdir = S.theta(ind)-phi; % main wave direction +# txtwdir = ['\theta_p=' num2pistr(Wdir,3)]; % convert to text string # -# Sf = spec2spec(S,'freq'); % frequency spectrum -# plotspec(Sf,1,lintype) +# plot([1 1]*S.theta(ind)-phi,[0 Dtheta(ind)],':'), hold on +# if LegendOn +# lh=legend(txtwdir,0); +# end +# plot(S.theta-phi,Dtheta,lintype) # -# subplot(212) +# fixthetalabels(thmin,thmax,'x',2) % fix xticklabel and xlabel for theta +# ylabel('D(\theta)') +# title('Spreading function') +# if ~ih, hold off, end +# %legend(lh) % refresh current legend +# elseif plotflag==4 % mesh +# mesh(freq,S.theta-phi,S.S) +# xlabel(xlbl_txt); +# fixthetalabels(thmin,thmax,'y',3) % fix yticklabel and ylabel for theta +# zlabel(zlbl_txt) +# title(ylbl3_txt) +# elseif plotflag==5 % mesh +# %h=polar([0 2*pi],[0 freq(end)]); +# %delete(h);hold on +# [X,Y]=meshgrid(S.theta-phi,freq); +# [X,Y]=polar2cart(X,Y); +# mesh(X,Y,S.S') +# % display the unit circle beneath the surface +# hold on, mesh(X,Y,zeros(size(S.S'))),hold off +# zlabel(zlbl_txt) +# title(ylbl3_txt) +# set(gca,'xticklabel','','yticklabel','') +# lighting phong +# %lighting gouraud +# %light +# elseif (plotflag==6) || (plotflag==7), +# theta = S.theta-phi; +# [c, h] = contour(freq,theta,S.S); %,Nlevel); % calculate levels +# fixthetalabels(thmin,thmax,'y',2) % fix yticklabel and ylabel for theta +# if plotflag==7, +# hold on +# [c,h] = contourf(freq,theta,S.S); %,Nlevel); % calculate levels +# %hold on +# end # -# Dtf = S.S; -# [Nt,Nf] = size(S.S); -# Sf = Sf.S(:).'; -# ind = find(Sf); +# title(ylbl3_txt) +# xlabel(xlbl_txt); +# if 0, +# [z_level] = clevels(c); % find contour levels +# % label the contour levels +# if txtFlag==1 +# textstart_x = 0.06; textstart_y=0.94; +# cltext1(z_level,textstart_x,textstart_y) % a local variant of cltext +# else +# cltext(z_level) +# end +# else +# colormap('jet') # -# if plotflag==3, %Directional distribution D(theta,freq)) -# Dtf(:,ind) = Dtf(:,ind)./Sf(ones(Nt,1),ind); -# end -# Dtheta = simpson(freq,Dtf,2); %Directional spreading, D(theta) -# Dtheta = Dtheta/simpson(S.theta,Dtheta); % make sure int D(theta)dtheta = 1 -# [y,ind] = max(Dtheta); -# Wdir = S.theta(ind)-phi; % main wave direction -# txtwdir = ['\theta_p=' num2pistr(Wdir,3)]; % convert to text string +# if plotflag==7, +# fcolorbar(c) +# else +# %clabel(c,h), +# hcb = colorbar; +# end +# grid on +# end +# else +# error('Unknown plot option') +# end +# otherwise, error('unknown spectral type') +# end # -# plot([1 1]*S.theta(ind)-phi,[0 Dtheta(ind)],':'), hold on -# if LegendOn -# lh=legend(txtwdir,0); -# end -# plot(S.theta-phi,Dtheta,lintype) +# if ~ih, hold off, end # -# fixthetalabels(thmin,thmax,'x',2) % fix xticklabel and xlabel for theta -# ylabel('D(\theta)') -# title('Spreading function') -# if ~ih, hold off, end -# %legend(lh) % refresh current legend -# elseif plotflag==4 % mesh -# mesh(freq,S.theta-phi,S.S) -# xlabel(xlbl_txt); -# fixthetalabels(thmin,thmax,'y',3) % fix yticklabel and ylabel for theta -# zlabel(zlbl_txt) -# title(ylbl3_txt) -# elseif plotflag==5 % mesh -# %h=polar([0 2*pi],[0 freq(end)]); -# %delete(h);hold on -# [X,Y]=meshgrid(S.theta-phi,freq); -# [X,Y]=polar2cart(X,Y); -# mesh(X,Y,S.S') -# % display the unit circle beneath the surface -# hold on, mesh(X,Y,zeros(size(S.S'))),hold off -# zlabel(zlbl_txt) -# title(ylbl3_txt) -# set(gca,'xticklabel','','yticklabel','') -# lighting phong -# %lighting gouraud -# %light -# elseif (plotflag==6) || (plotflag==7), -# theta = S.theta-phi; -# [c, h] = contour(freq,theta,S.S); %,Nlevel); % calculate levels -# fixthetalabels(thmin,thmax,'y',2) % fix yticklabel and ylabel for theta -# if plotflag==7, -# hold on -# [c,h] = contourf(freq,theta,S.S); %,Nlevel); % calculate levels -# %hold on -# end +# # The following two commands install point-and-click editing of +# # all the text objects (title, xlabel, ylabel) of the current figure: # -# title(ylbl3_txt) -# xlabel(xlbl_txt); -# if 0, -# [z_level] = clevels(c); % find contour levels -# % label the contour levels -# if txtFlag==1 -# textstart_x = 0.06; textstart_y=0.94; -# cltext1(z_level,textstart_x,textstart_y) % a local variant of cltext -# else -# cltext(z_level) -# end -# else -# colormap('jet') -# -# if plotflag==7, -# fcolorbar(c) -# else -# %clabel(c,h), -# hcb = colorbar; -# end -# grid on -# end -# else -# error('Unknown plot option') -# end -# otherwise, error('unknown spectral type') -# end -# -# if ~ih, hold off, end -# -# % The following two commands install point-and-click editing of -# % all the text objects (title, xlabel, ylabel) of the current figure: -# -# %set(findall(gcf,'type','text'),'buttondownfcn','edtext') -# %set(gcf,'windowbuttondownfcn','edtext(''hide'')') -# -# return -# - - -# -# function fixthetalabels(thmin,thmax,xy,dim) -# %FIXTHETALABELS pretty prints the ticklabels and x or y labels for theta -# % -# % CALL fixthetalabels(thmin,thmax,xy,dim) -# % -# % thmin, thmax = minimum and maximum value for theta (wave direction) -# % xy = 'x' if theta is plotted on the x-axis -# % 'y' if theta is plotted on the y-axis -# % dim = specifies the dimension of the plot (ie number of axes shown 2 or 3) -# % If abs(thmax-thmin)<3*pi it is assumed that theta is given in radians -# % otherwise degrees -# -# ind = [('x' == xy) ('y' == xy) ]; -# yx = 'yx'; -# yx = yx(ind); -# if nargin<4||isempty(dim), -# dim=2; -# end -# %drawnow -# %pause -# -# if abs(thmax-thmin)<3*pi, %Radians given. Want xticks given as fractions of pi -# %Trick to update the axis -# if xy=='x' -# if dim<3, -# axis([thmin,thmax 0 inf ]) -# else -# axis([thmin,thmax 0 inf 0 inf]) -# end -# else -# if dim<3, -# axis([0 inf thmin,thmax ]) -# else -# axis([0 inf thmin,thmax 0 inf]) -# end -# end -# -# set(gca,[xy 'tick'],pi*(thmin/pi:0.25:thmax/pi)); -# set(gca,[xy 'ticklabel'],[]); -# x = get(gca,[xy 'tick']); -# y = get(gca,[yx 'tick']); -# y1 = y(1); -# dy = y(2)-y1; -# yN = y(end)+dy; -# ylim = [y1 yN]; -# dy1 = diff(ylim)/40; -# %ylim=get(gca,[yx 'lim'])%,ylim=ylim(2); -# -# if xy=='x' -# for j=1:length(x) -# xtxt = num2pistr(x(j)); -# figtext(x(j),y1-dy1,xtxt,'data','data','center','top'); -# end -# % ax = [thmin thmax 0 inf]; -# ax = [thmin thmax ylim]; -# if dim<3, -# figtext(mean(x),y1-7*dy1,'Wave directions (rad)','data','data','center','top') -# else -# ax = [ax 0 inf]; -# xlabel('Wave directions (rad)') -# end -# else -# %ax = [0 inf thmin thmax]; -# ax = [ylim thmin thmax]; -# -# if dim<3, -# for j=1:length(x) -# xtxt = num2pistr(x(j)); -# figtext(y1-dy1/2,x(j),xtxt,'data','data','right'); -# end -# set(gca,'DefaultTextRotation',90) -# %ylabel('Wave directions (rad)') -# figtext(y1-3*dy1,mean(x),'Wave directions (rad)','data','data','center','bottom') -# set(gca,'DefaultTextRotation',0) -# else -# for j=1:length(x) -# xtxt = num2pistr(x(j)); -# figtext(y1-3*dy1,x(j),xtxt,'data','data','right'); -# end -# ax = [ax 0 inf]; -# ylabel('Wave directions (rad)') -# end -# end -# %xtxt = num2pistr(x(j)); -# %for j=2:length(x) -# % xtxt = strvcat(xtxt,num2pistr(x(j))); -# %end -# %set(gca,[xy 'ticklabel'],xtxt) -# else % Degrees given -# set(gca,[xy 'tick'],thmin:45:thmax) -# if xy=='x' -# ax=[thmin thmax 0 inf]; -# if dim>=3, ax=[ax 0 inf]; end -# xlabel('Wave directions (deg)') -# else -# ax=[0 inf thmin thmax ]; -# if dim>=3, ax=[ax 0 inf]; end -# ylabel('Wave directions (deg)') -# end -# end -# axis(ax) -# return +# #set(findall(gcf,'type','text'),'buttondownfcn','edtext') +# #set(gcf,'windowbuttondownfcn','edtext(''hide'')') # +# return class SpecData1D(PlotData): @@ -684,6 +573,43 @@ class SpecData1D(PlotData): self.setlabels() + def _get_default_dt_and_rate(self, dt): + dt_old = self.sampling_period() + if dt is None: + return dt_old, 1 + rate = max(round(dt_old * 1. / dt), 1.) + return dt, rate + + def _check_dt(self, dt): + freq = self.args + checkdt = 1.2 * min(diff(freq)) / 2. / pi + if self.freqtype in 'f': + checkdt *= 2 * pi + if (checkdt < 2. ** -16 / dt): + print('Step dt = %g in computation of the density is ' + + 'too small.' % dt) + print('The computed covariance (by FFT(2^K)) may differ from the') + print('theoretical. Solution:') + raise ValueError('use larger dt or sparser grid for spectrum.') + + def _check_cov_matrix(self, acfmat, nt, dt): + eps0 = 0.0001 + if nt + 1 >= 5: + cc2 = acfmat[0, 0] - acfmat[4, 0] * (acfmat[4, 0] / acfmat[0, 0]) + if (cc2 < eps0): + warnings.warn('Step dt = %g in computation of the density ' + + 'is too small.' % dt) + cc1 = acfmat[0, 0] - acfmat[1, 0] * (acfmat[1, 0] / acfmat[0, 0]) + if (cc1 < eps0): + warnings.warn('Step dt = %g is small, and may cause numerical ' + + 'inaccuracies.' % dt) + + @property + def lagtype(self): + if self.freqtype in 'k': # options are 'f' and 'w' and 'k' + return 'x' + return 't' + def tocov_matrix(self, nr=0, nt=None, dt=None): ''' Computes covariance function and its derivatives, alternative version @@ -725,59 +651,30 @@ class SpecData1D(PlotData): objects ''' - ftype = self.freqtype # %options are 'f' and 'w' and 'k' + + dt, rate = self._get_default_dt_and_rate(dt) + self._check_dt(dt) + freq = self.args n_f = len(freq) - dt_old = self.sampling_period() - if dt is None: - dt = dt_old - rate = 1 - else: - rate = max(round(dt_old * 1. / dt), 1.) - if nt is None: nt = rate * (n_f - 1) else: # %check if Nt is ok nt = minimum(nt, rate * (n_f - 1)) - checkdt = 1.2 * min(diff(freq)) / 2. / pi - if ftype in 'k': - lagtype = 'x' - else: - lagtype = 't' - if ftype in 'f': - checkdt = checkdt * 2 * pi - msg1 = 'Step dt = %g in computation of the density is too small.' % dt - msg2 = 'Step dt = %g is small, and may cause numerical inaccuracies.' % dt - - if (checkdt < 2. ** -16 / dt): - print(msg1) - print('The computed covariance (by FFT(2^K)) may differ from the') - print('theoretical. Solution:') - raise ValueError('use larger dt or sparser grid for spectrum.') - - # Calculating covariances - #~~~~~~~~~~~~~~~~~~~~~~~~ spec = self.copy() spec.resample(dt) acf = spec.tocovdata(nr, nt, rate=1) acfmat = zeros((nt + 1, nr + 1), dtype=float) acfmat[:, 0] = acf.data[0:nt + 1] - fieldname = 'R' + lagtype * nr + fieldname = 'R' + self.lagtype * nr for i in range(1, nr + 1): fname = fieldname[:i + 1] r_i = getattr(acf, fname) acfmat[:, i] = r_i[0:nt + 1] - eps0 = 0.0001 - if nt + 1 >= 5: - cc2 = acfmat[0, 0] - acfmat[4, 0] * (acfmat[4, 0] / acfmat[0, 0]) - if (cc2 < eps0): - warnings.warn(msg1) - cc1 = acfmat[0, 0] - acfmat[1, 0] * (acfmat[1, 0] / acfmat[0, 0]) - if (cc1 < eps0): - warnings.warn(msg2) + self._check_cov_matrix(acfmat, nt, dt) return acfmat def tocovdata(self, nr=0, nt=None, rate=None): @@ -1657,8 +1554,8 @@ class SpecData1D(PlotData): # ftmp = cov2mmtpdfexe(R,dt,u,defnr,Nstart,hg,options) # err = repmat(nan,size(ftmp)) # else - [ftmp, err, terr, options] = self._cov2mmtpdf( - R, dt, u, defnr, Nstart, hg, options) + ftmp, err, terr, options = self._cov2mmtpdf(R, dt, u, defnr, Nstart, + hg, options) # end note = '' @@ -1850,8 +1747,8 @@ class SpecData1D(PlotData): (indI(3)=Nt+1); for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn)) (indI(4)=Nt+2); for i\in (indI(4)+1,indI(5)], Y(i)<0 (deriv. X'(ts)) ''' - R0, R1, R2, R3, R4, R5 = R[:, :5].T - + R0, R1, R2, R3, R4 = R[:, :5].T + covinput = self._covinput_mmt_pdf Ntime = len(R0) Nx0 = max(1, len(hg)) Nx1 = Nx0 @@ -1876,14 +1773,15 @@ class SpecData1D(PlotData): if def_nr <= 1: # just plain Mm Nx = Nx1 * (Nx1 - 1) / 2 IJ = (Nx1 + isOdd) / 2 - if (hg[0] + hg[Nx1 - 1] == 0 and (hg[IJ - 1] == 0 or hg[IJ - 1] + hg[IJ] == 0)): + if (hg[0] + hg[Nx1 - 1] == 0 and (hg[IJ - 1] == 0 or + hg[IJ - 1] + hg[IJ] == 0)): symmetry = 0 print(' Integration region symmetric') # May save Nx1-isOdd integrations in each time step # This is not implemented yet. # Nx = Nx1*(Nx1-1)/2-Nx1+isOdd - - # CC = normalizing constant = 1/ expected number of zero-up-crossings of X' + # normalizing constant: + # CC = 1/ expected number of zero-up-crossings of X' # CC = 2*pi*sqrt(-R2[0]/R4[0]) # XcScale = log(CC) XcScale = log(2 * pi * sqrt(-R2[0] / R4[0])) @@ -1902,16 +1800,15 @@ class SpecData1D(PlotData): Nc = 5 NI = 5 Nd = 3 - # CC= normalizing constant= 1/ expected number of u-up-crossings of X + # CC = 1/ expected number of u-up-crossings of X # CC = 2*pi*sqrt(-R0(1)/R2(1))*exp(0.5D0*u*u/R0(1)) XcScale = log(2 * pi * sqrt(-R0[0] / R2[0])) + 0.5 * u * u / R0[0] options['xcscale'] = XcScale - opt0 = struct2cell(options) - # opt0 = opt0(1:10) - # seed = [] - # opt0 = {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT1} - +# opt0 = [options[n] for n in ('SCIS', 'XcScale', 'ABSEPS', 'RELEPS', +# 'COVEPS', 'MAXPTS', 'MINPTS', 'seed', +# 'NIT1')] + rind = Rind(**options) if (Nx > 1): # (M,m) or (M,m)v distribution wanted if ((def_nr == 0 or def_nr == 2)): @@ -1982,10 +1879,10 @@ class SpecData1D(PlotData): indI[2] = Nt + 1 indI[3] = Ntd # positive wave period - BIG[:Ntdc, :Ntdc] = covinput( - BIG[:Ntdc, :Ntdc], R0, R1, R2, R3, R4, Ntd, 0) + # self._covinput_mmt_pdf(BIG, R, tn, ts, tnold) + BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], R, Ntd, 0) [fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], - a_lo, a_up, indI, xc, Nt, *opt0) + a_lo, a_up, indI, xc, Nt) # fxind = CC*rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,Nt,NIT1, # speed1,indI,a_lo,a_up) if (Nx < 2): @@ -2001,12 +1898,9 @@ class SpecData1D(PlotData): if def_nr in [-2, -1, 0]: for i in range(1, Nx1): J = IJ + i - pdf[:i, i, 0] = pdf[:i, i, 0] + \ - fxind[IJ:J].T * dt # *CC - err[:i, i, 0] = err[:i, i, 0] + \ - (err0[IJ + 1:J].T * dt) ** 2 - terr[:i, i, 0] = terr[ - :i, i, 0] + (terr0[IJ:J].T * dt) + pdf[:i, i, 0] += fxind[IJ:J].T * dt # *CC + err[:i, i, 0] += (err0[IJ + 1:J].T * dt) ** 2 + terr[:i, i, 0] += (terr0[IJ:J].T * dt) IJ = J elif def_nr == 1: # joint density of (M,m,TMm) for i in range(1, Nx1): @@ -2020,24 +1914,18 @@ class SpecData1D(PlotData): elif def_nr == 2: for i in range(1, Nx1): J = IJ + Nx1 - pdf[1:Nx1, i, 0] = pdf[1:Nx1, i, 0] + \ - fxind[IJ:J].T * dt # %*CC - err[1:Nx1, i, 0] = err[1:Nx1, i, 0] + \ - (err0[IJ:J].T * dt) ** 2 - terr[1:Nx1, i, 0] = terr[ - 1:Nx1, i, 0] + (terr0[IJ:J].T * dt) + pdf[1:Nx1, i, 0] += fxind[IJ:J].T * dt # %*CC + err[1:Nx1, i, 0] += (err0[IJ:J].T * dt) ** 2 + terr[1:Nx1, i, 0] += (terr0[IJ:J].T * dt) IJ = J # end %do elif def_nr == 3: # joint density of level v separated (M,m,TMm)v for i in range(1, Nx1): J = IJ + Nx1 - pdf[1:Nx1, i, Ntd] = pdf[ - 1:Nx1, i, Ntd] + fxind[IJ:J].T # %*CC - err[1:Nx1, i, Ntd] = err[ - 1:Nx1, i, Ntd] + (err0[IJ:J].T) ** 2 - terr[1:Nx1, i, Ntd] = terr[ - 1:Nx1, i, Ntd] + (terr0[IJ:J].T) + pdf[1:Nx1, i, Ntd] += fxind[IJ:J].T # %*CC + err[1:Nx1, i, Ntd] += (err0[IJ:J].T) ** 2 + terr[1:Nx1, i, Ntd] += (terr0[IJ:J].T) IJ = J # end %do # end % SELECT @@ -2067,10 +1955,9 @@ class SpecData1D(PlotData): for ts in range(1, tn - 1): # = 2:tn-1: # positive wave period BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], - R0, R1, R2, R3, R4, - tn, ts, tnold) - [fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], - a_lo, a_up, indI, xc, Nt, *opt0) + R, tn, ts, tnold) + fxind, err0, terr0 = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], + a_lo, a_up, indI, xc, Nt) # tnold = tn if def_nr in [3, 4]: @@ -2087,12 +1974,9 @@ class SpecData1D(PlotData): IJ = 0 for i in range(1, Nx1): J = IJ + Nx1 - pdf[1:Nx1, i, ts] = pdf[ - 1:Nx1, i, ts] + fxind[IJ:J].T * dt - err[1:Nx1, i, ts] = err[ - 1:Nx1, i, ts] + (err0[IJ:J].T * dt) ** 2 - terr[1:Nx1, i, ts] = terr[ - 1:Nx1, i, ts] + (terr0[IJ:J].T * dt) + pdf[1:Nx1, i, ts] += fxind[IJ:J].T * dt + err[1:Nx1, i, ts] += (err0[IJ:J].T * dt) ** 2 + terr[1:Nx1, i, ts] += (terr0[IJ:J].T * dt) IJ = J # end %do # end @@ -2111,12 +1995,9 @@ class SpecData1D(PlotData): for i in range(1, Nx1): # = 2:Nx1 J = IJ + Nx1 # %*CC - pdf[1:Nx1, i, tn - ts] = pdf[1:Nx1, - i, tn - ts] + fxind[IJ:J].T * dt - err[1:Nx1, i, tn - ts] = err[1:Nx1, i, - tn - ts] + (err0[IJ:J].T * dt) ** 2 - terr[ - 1:Nx1, i, tn - ts] = terr[1:Nx1, i, tn - ts] + (terr0[IJ:J].T * dt) + pdf[1:Nx1, i, tn - ts] += fxind[IJ:J].T * dt + err[1:Nx1, i, tn - ts] += (err0[IJ:J].T * dt) ** 2 + terr[1:Nx1, i, tn - ts] += (terr0[IJ:J].T * dt) IJ = J # end %do # end @@ -2129,10 +2010,9 @@ class SpecData1D(PlotData): # linear. # positive wave period BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], - R0, R1, R2, R3, R4, - tn, ts, tnold) - [fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], - a_lo, a_up, indI, xc, Nt, *opt0) + R, tn, ts, tnold) + fxind, err0, terr0 = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], + a_lo, a_up, indI, xc, Nt) #[fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex,a_lo,a_up,indI, xc,Nt,opt0{:}) # tnold = tn @@ -2144,9 +2024,9 @@ class SpecData1D(PlotData): err[0, ts, tn] = err0[0] ** 2 err[0, ts, tn] = terr0(1) if (ts < tn - ts): # %THEN - pdf[0, tn - ts, tn] = fxind[0] # *CC - err[0, tn - ts, tn] = err0[0] ** 2 - terr[0, tn - ts, tn] = terr0[0] + pdf[0, tn - ts, tn] = fxind[0] # *CC + err[0, tn - ts, tn] = err0[0] ** 2 + terr[0, tn - ts, tn] = terr0[0] # end # GOTO 350 else: @@ -2156,43 +2036,32 @@ class SpecData1D(PlotData): # Max to the crossing of level u (M,m,TMd). for i in range(1, Nx1): J = IJ + Nx1 - pdf[1:Nx1, i, ts] = pdf[ - 1:Nx1, i, ts] + fxind[IJ:J] * dt # %*CC - err[1:Nx1, i, ts] = err[ - 1:Nx1, i, ts] + (err0[IJ:J] * dt) ** 2 - terr[1:Nx1, i, ts] = terr[ - 1:Nx1, i, ts] + (terr0[IJ:J] * dt) + pdf[1:Nx1, i, ts] += fxind[IJ:J] * dt # %*CC + err[1:Nx1, i, ts] += (err0[IJ:J] * dt) ** 2 + terr[1:Nx1, i, ts] += (terr0[IJ:J] * dt) if (ts < tn - ts): # exploiting the symmetry # %*CC - pdf[i, 1:Nx1, tn - ts] = pdf[i, - 1:Nx1, tn - ts] + fxind[IJ:J] * dt - err[i, 1:Nx1, tn - ts] = err[i, 1:Nx1, - tn - ts] + (err0[IJ:J] * dt) ** 2 - terr[ - i, 1:Nx1, tn - ts] = terr[i, 1:Nx1, tn - ts] + (terr0[IJ:J] * dt) + pdf[i, 1:Nx1, tn - ts] += fxind[IJ:J] * dt + err[i, 1:Nx1, tn - ts] += (err0[IJ:J] * dt) ** 2 + terr[i, 1:Nx1, tn - ts] += (terr0[IJ:J] * dt) # end IJ = J - # end %do + # end %do elif def_nr == 5: # 5, gives level u separated Max2min and wave period # from the crossing of level u to min (M,m,Tdm). for i in range(1, Nx1): # = 2:Nx1, J = IJ + Nx1 - pdf[1:Nx1, i, tn - ts] = pdf[1:Nx1, - i, tn - ts] + fxind[IJ:J] * dt - err[1:Nx1, i, tn - ts] = err[1:Nx1, i, - tn - ts] + (err0[IJ:J] * dt) ** 2 + pdf[1:Nx1, i, tn - ts] += fxind[IJ:J] * dt + err[1:Nx1, i, tn - ts] += (err0[IJ:J] * dt) ** 2 terr[ - 1:Nx1, i, tn - ts] = terr[1:Nx1, i, tn - ts] + (terr0[IJ:J] * dt) + 1:Nx1, i, tn - ts] += (terr0[IJ:J] * dt) if (ts < tn - ts + 1): # exploiting the symmetry - pdf[i, 1:Nx1, ts] = pdf[ - i, 1:Nx1, ts] + fxind[IJ:J] * dt - err[i, 1:Nx1, ts] = err[ - i, 1:Nx1, ts] + (err0[IJ:J] * dt) ** 2 - terr[i, 1:Nx1, ts] = terr[ - i, 1:Nx1, ts] + (terr0[IJ:J] * dt) + pdf[i, 1:Nx1, ts] += fxind[IJ:J] * dt + err[i, 1:Nx1, ts] += (err0[IJ:J] * dt) ** 2 + terr[i, 1:Nx1, ts] += (terr0[IJ:J] * dt) # end %ENDIF IJ = J # end %do @@ -2263,7 +2132,7 @@ class SpecData1D(PlotData): Cov(X''(t),X(s)) = r''(s-t) = r''(|s-t|) Cov(X''(t),X''(s)) = r''''(s-t) = r''''(|s-t|) """ - R0, R1, R2, R3, R4 = R.T + R0, R1, R2, R3, R4 = R[:, :5].T if (ts > 1): shft = 1 N = tn + 5 + shft @@ -4006,17 +3875,17 @@ class SpecData2D(PlotData): # ftype = self.freqtype freq = self.args[1] theta = linspace(-pi, pi, ntOld) - [F, T] = meshgrid(freq, theta) + # [F, T] = meshgrid(freq, theta) dtheta = self.theta[1] - self.theta[0] self.theta[nt] = self.theta[nt - 1] + dtheta self.data[nt, :] = self.data[0, :] - self.data = interp2(freq, - np.vstack([self.theta[0] - dtheta, - self.theta]), - np.vstack([self.data[nt, :], - self.data]), F, T, - method) + self.data = interp2d(freq, + np.vstack([self.theta[0] - dtheta, + self.theta]), + np.vstack([self.data[nt, :], + self.data]), + kind=method)(freq, theta) self.args[0] = theta elif stype == 'k2d': @@ -4274,33 +4143,16 @@ class SpecData2D(PlotData): self.labels.zlab = labels[2] -def _test_specdata(): - import wafo.spectrum.models as sm - Sj = sm.Jonswap() - S = Sj.tospecdata() - me, va, sk, ku = S.stats_nl(moments='mvsk') - - def main(): import matplotlib matplotlib.interactive(True) from wafo.spectrum import models as sm - w = linspace(0, 3, 100) Sj = sm.Jonswap() S = Sj.tospecdata() - f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7) - f.err - f.plot() - f.show() - # pdfplot(f) - # hold on, - # plot(f.x{:}, f.f+f.err,'r',f.x{:}, f.f-f.err) estimated error bounds - # hold off - # S = SpecData1D(Sj(w),w) R = S.tocovdata(nr=1) - S1 = S.copy() + Si = R.tospecdata() ns = 5000 dt = .2 diff --git a/wafo/spectrum/models.py b/wafo/spectrum/models.py index 03d2bc0..45619aa 100644 --- a/wafo/spectrum/models.py +++ b/wafo/spectrum/models.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python """ Models module ------------- @@ -26,17 +27,7 @@ Spreading - Directional spreading function. """ -# Name: models -# Purpose: Interface to various spectrum models -# -# Author: pab -# -# Created: 29.08.2008 -# Copyright: (c) pab 2008 -# Licence: - -#!/usr/bin/env python -from __future__ import division +from __future__ import absolute_import, division import warnings from scipy.interpolate import interp1d @@ -50,17 +41,20 @@ from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, cos, abs, sinh, isfinite, mod, expm1, tanh, cosh, finfo, ones, ones_like, isnan, zeros_like, flatnonzero, sinc, hstack, vstack, real, flipud, clip) -from wafo.wave_theory.dispersion_relation import w2k, k2w # @UnusedImport -from wafo.spectrum import SpecData1D, SpecData2D -sech = lambda x: 1.0 / cosh(x) - -eps = finfo(float).eps +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 @@ -409,7 +403,7 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, # The following formulas are from Lewis and Allos 1990: zeta = g * fetch / (u10 ** 2) # dimensionless fetch, Table 1 - #zeta = min(zeta, 2.414655013429281e+004) + # zeta = min(zeta, 2.414655013429281e+004) if method.startswith('h'): if method[-1] == '3': # Hasselman et.al (1973) A = 0.076 * zeta ** (-0.22) @@ -543,7 +537,7 @@ class Jonswap(ModelSpectrum): self.method = method self.wnc = wnc - if self.gamma == None or not isfinite(self.gamma) or self.gamma < 1: + if self.gamma is None or not isfinite(self.gamma) or self.gamma < 1: self.gamma = jonswap_peakfact(Hm0, Tp) self._preCalculateAg() @@ -580,7 +574,7 @@ class Jonswap(ModelSpectrum): if self.gamma == 1: self.Ag = 1.0 self.method = 'parametric' - elif self.Ag != None: + elif self.Ag is not None: self.method = 'custom' if self.Ag <= 0: raise ValueError('Ag must be larger than 0!') @@ -615,7 +609,7 @@ class Jonswap(ModelSpectrum): # elseif N == 5 && M == 4, # options.Ag = (1+1.0*log(gammai).**1.16)/gammai # options.Ag = (1-0.287*log(gammai)) - ### options.normalizeMethod = 'Three' + # options.normalizeMethod = 'Three' # elseif N == 4 && M == 4, # options.Ag = (1+1.1*log(gammai).**1.19)/gammai else: @@ -624,7 +618,7 @@ class Jonswap(ModelSpectrum): 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') + 'sigmaA!=0.07 or sigmaB!=0.09') def peak_e_factor(self, wn): ''' PEAKENHANCEMENTFACTOR @@ -790,9 +784,9 @@ class Tmaspec(Jonswap): self.type = 'TMA' def phi(self, w, h=None, g=None): - if h == None: + if h is None: h = self.h - if g == None: + if g is None: g = self.g return phi1(w, h, g) @@ -1005,8 +999,8 @@ class Torsethaugen(ModelSpectrum): 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)) + # 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) ** @@ -1014,7 +1008,7 @@ class Torsethaugen(ModelSpectrum): else: Tpw = inf - #Tpw = max(Tpw,2.5) + # Tpw = max(Tpw,2.5) gammaw = 1 if monitor: if Rpw > 0.1: @@ -1095,14 +1089,14 @@ class McCormick(Bretschneider): self.type = 'McCormick' self.Hm0 = Hm0 self.Tp = Tp - if Tz == None: + if Tz is None: Tz = 0.8143 * Tp self.Tz = Tz if chk_seastate: self.chk_seastate() - if M == None and self.Hm0 > 0: + if M is None and self.Hm0 > 0: self._TpdTz = Tp / Tz M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5) self.M = M @@ -1410,7 +1404,7 @@ class Spreading(object): 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 + 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 ] @@ -1510,7 +1504,7 @@ class Spreading(object): if not self.method[0] in methods: raise ValueError('Unknown method') self.method = methods[self.method[0]] - elif self.method == None: + elif self.method is None: pass else: if method < 0 or 3 < method: @@ -1561,7 +1555,7 @@ class Spreading(object): '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 sqrt(eps) * abs(A)) * (abs(da) > sqrt(eps))) + (abs(da) > sqrt(_EPS) * abs(A)) * (abs(da) > sqrt(_EPS))) if ix.size == 0: if any(A > pi): raise ValueError( @@ -1837,8 +1831,10 @@ class Spreading(object): ''' 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))) - fun0 = lambda x: sp.ive(1, x) / sp.ive(0, x) funK = interp1d(fun0(K0), K0) K0 = funK(r1.ravel()) k1 = flatnonzero(isnan(K0)) @@ -1848,15 +1844,14 @@ class Spreading(object): ix0 = flatnonzero(r1 != 0.0) K = zeros_like(r1) - fun = lambda x: fun0(x) - r1[ix] for ix in ix0: - K[ix] = optimize.fsolve(fun, K0[ix]) + 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))) + B0 = hstack((linspace(_EPS, 5, 513), linspace(5.0001, 100))) funB = interp1d(self._r1ofsech2(B0), B0) B0 = funB(r1.ravel()) @@ -1867,7 +1862,9 @@ class Spreading(object): ix0 = flatnonzero(r1 != 0.0) B = zeros_like(r1) - fun = lambda x: 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix] + + def fun(x): + return 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix] for ix in ix0: B[ix] = abs(optimize.fsolve(fun, B0[ix])) return B @@ -1890,7 +1887,7 @@ class Spreading(object): S : ndarray spread parameter of COS2S functions ''' - if self.method == None: + if self.method is None: # no frequency dependent spreading, # but possible frequency dependent direction s = atleast_1d(self.s_a) @@ -1923,12 +1920,12 @@ class Spreading(object): # 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 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 + # Convert to S-paramater in COS-2S distribution s = r1 / (1. - r1) else: s[k] = 0.0 @@ -1994,7 +1991,7 @@ class Spreading(object): theta = np.linspace(-pi, pi, nt) else: L = abs(theta[-1] - theta[0]) - if abs(L - pi) > eps: + if abs(L - pi) > _EPS: raise ValueError('theta must cover all angles -pi -> pi') nt = len(theta) @@ -2015,7 +2012,7 @@ class Spreading(object): Snew.h = specdata.h Snew.phi = phi0 Snew.norm = specdata.norm - #Snew.note = specdata.note + ', spreading: %s' % self.type + # Snew.note = specdata.note + ', spreading: %s' % self.type return Snew @@ -2036,11 +2033,9 @@ def _test_some_spectra(): plb.show() plb.close('all') - #import pylab as plb - #w = plb.linspace(0,3) w, th = plb.ogrid[0:4, 0:6] k, k2 = w2k(w, th) - #k1, k12 = w2k(w, th, h=20) + plb.plot(w, k, w, k2) plb.show() @@ -2080,8 +2075,8 @@ def _test_spreading(): pi = plb.pi w = plb.linspace(0, 3, 257) theta = plb.linspace(-pi, pi, 129) - theta0 = lambda w: w * plb.pi / 6.0 - D2 = Spreading('cos2s', theta0=theta0) + + D2 = Spreading('cos2s', theta0=lambda w: w * plb.pi / 6.0) d1 = D2(theta, w)[0] _t = plb.contour(d1.squeeze()) diff --git a/wafo/spectrum/tests/__init__.py b/wafo/spectrum/tests/__init__.py index 8b13789..e69de29 100644 --- a/wafo/spectrum/tests/__init__.py +++ b/wafo/spectrum/tests/__init__.py @@ -1 +0,0 @@ - diff --git a/wafo/spectrum/tests/test_specdata1d.py b/wafo/spectrum/tests/test_specdata1d.py index 3a7e900..645c1fa 100644 --- a/wafo/spectrum/tests/test_specdata1d.py +++ b/wafo/spectrum/tests/test_specdata1d.py @@ -3,7 +3,8 @@ import wafo.transform.models as wtm import wafo.objects as wo from wafo.spectrum import SpecData1D import numpy as np -from numpy.testing import assert_array_almost_equal +from numpy import NAN +from numpy.testing import assert_array_almost_equal, assert_array_equal import unittest @@ -12,213 +13,184 @@ def slow(f): return f -class TestSpectrum(unittest.TestCase): +class TestSpectrumHs7(unittest.TestCase): + def setUp(self): + self.Sj = sm.Jonswap(Hm0=7.0, Tp=11) + self.S = self.Sj.tospecdata() - @slow def test_tocovmatrix(self): - Sj = sm.Jonswap() - S = Sj.tospecdata() - acfmat = S.tocov_matrix(nr=3, nt=256, dt=0.1) + acfmat = self.S.tocov_matrix(nr=3, nt=256, dt=0.1) vals = acfmat[:2, :] true_vals = np.array([[3.06073383, 0.0000000, -1.67748256, 0.], [3.05235423, -0.1674357, -1.66811444, 0.18693242]]) - self.assertTrue((np.abs(vals - true_vals) < 1e-7).all()) - - -def test_tocovdata(): - Sj = sm.Jonswap() - S = Sj.tospecdata() - Nt = len(S.data) - 1 - acf = S.tocovdata(nr=0, nt=Nt) - vals = acf.data[:5] - - true_vals = np.array( - [3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042]) - assert((np.abs(vals - true_vals) < 1e-6).all()) - - -def test_to_t_pdf(): - Sj = sm.Jonswap() - S = Sj.tospecdata() - f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100) - vals = ['%2.3f' % val for val in f.data[:10]] - truevals = ['0.000', '0.014', '0.027', '0.040', - '0.050', '0.059', '0.067', '0.073', '0.077', '0.082'] - for t, v in zip(truevals, vals): - assert(t == v) - - # estimated error bounds - vals = ['%2.4f' % val for val in f.err[:10]] - truevals = ['0.0000', '0.0003', '0.0003', '0.0004', - '0.0006', '0.0008', '0.0016', '0.0019', '0.0020', '0.0021'] - for t, v in zip(truevals, vals): - assert(t == v) - - -@slow -def test_sim(): - Sj = sm.Jonswap() - S = Sj.tospecdata() - #ns = 100 - #dt = .2 - #x1 = S.sim(ns, dt=dt) - - import scipy.stats as st - x2 = S.sim(20000, 20) - truth1 = [0, np.sqrt(S.moment(1)[0]), 0., 0.] - funs = [np.mean, np.std, st.skew, st.kurtosis] - for fun, trueval in zip(funs, truth1): - res = fun(x2[:, 1::], axis=0) - m = res.mean() - sa = res.std() - #trueval, m, sa - assert(np.abs(m - trueval) < sa) - - -@slow -def test_sim_nl(): - - Sj = sm.Jonswap() - S = Sj.tospecdata() -# ns = 100 -# dt = .2 -# x1 = S.sim_nl(ns, dt=dt) - import scipy.stats as st - x2, _x1 = S.sim_nl(ns=20000, cases=40) - truth1 = [0, np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') - truth1[-1] = truth1[-1] - 3 - - # truth1 - #[0, 1.7495200310090633, 0.18673120577479801, 0.061988521262417606] - - funs = [np.mean, np.std, st.skew, st.kurtosis] - for fun, trueval in zip(funs, truth1): - res = fun(x2.data, axis=0) - m = res.mean() - sa = res.std() - # trueval, m, sa - assert(np.abs(m - trueval) < 2 * sa) - - -def test_stats_nl(): - - Hs = 7. - Sj = sm.Jonswap(Hm0=Hs, Tp=11) - S = Sj.tospecdata() - me, va, sk, ku = S.stats_nl(moments='mvsk') - assert(me == 0.0) - assert_array_almost_equal(va, 3.0608203389019537) - assert_array_almost_equal(sk, 0.18673120577479801) - assert_array_almost_equal(ku, 3.0619885212624176) - - -def test_testgaussian(): - - Hs = 7 - Sj = sm.Jonswap(Hm0=Hs) - S0 = Sj.tospecdata() - # ns =100; dt = .2 - # x1 = S0.sim(ns, dt=dt) - - S = S0.copy() - me, _va, sk, ku = S.stats_nl(moments='mvsk') - S.tr = wtm.TrHermite( - mean=me, sigma=Hs / 4, skew=sk, kurt=ku, ysigma=Hs / 4) - ys = wo.mat2timeseries(S.sim(ns=2 ** 13)) - g0, _gemp = ys.trdata() - t0 = g0.dist2gauss() - t1 = S0.testgaussian(ns=2 ** 13, test0=None, cases=50) - assert(sum(t1 > t0) < 5) - - -def test_moment(): - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() # Make spectrum ob - vals, txt = S.moment() - true_vals = [1.5614600345079888, 0.95567089481941048] - true_txt = ['m0', 'm0tt'] - for tv, v in zip(true_vals, vals): - assert_array_almost_equal(tv, v) - for tv, v in zip(true_txt, txt): - assert(tv==v) - - -def test_nyquist_freq(): - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() # Make spectrum ob - assert(S.nyquist_freq() == 3.0) - - -def test_sampling_period(): - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() # Make spectrum ob - assert(S.sampling_period() == 1.0471975511965976) - - -def test_normalize(): - Sj = sm.Jonswap(Hm0=5) - S = Sj.tospecdata() # Make spectrum ob - S.moment(2) - ([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt']) - vals, _txt = S.moment(2) - true_vals = [1.5614600345079888, 0.95567089481941048] - for tv, v in zip(true_vals, vals): - assert_array_almost_equal(tv, v) - - Sn = S.copy() - Sn.normalize() - - # Now the moments should be one - new_vals, _txt = Sn.moment(2) - for v in new_vals: - assert(np.abs(v - 1.0) < 1e-7) - - -def test_characteristic(): - ''' - >>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap(Hm0=5) - >>> S = Sj.tospecdata() #Make spectrum ob - >>> S.characteristic(1) - (array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) - - >>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers - >>> ch; R; txt - array([ 8.59007646, 8.03139757, 5.62484314]) - array([[ 0.03040216, 0.02834263, nan], - [ 0.02834263, 0.0274645 , nan], - [ nan, nan, 0.01500249]]) - ['Tm01', 'Tm02', 'Tm24'] - - >>> S.characteristic('Ss') # fact a string - (array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) - - >>> S.characteristic(['Hm0','Tm02']) # fact a list of strings - (array([ 4.99833578, 8.03139757]), array([[ 0.05292989, 0.02511371], - [ 0.02511371, 0.0274645 ]]), ['Hm0', 'Tm02']) - ''' - + assert_array_almost_equal(vals, true_vals) + + def test_tocovdata(self): + + Nt = len(self.S.data) - 1 + acf = self.S.tocovdata(nr=0, nt=Nt) + vals = acf.data[:5] + + true_vals = np.array( + [3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042]) + assert_array_almost_equal(vals, true_vals) + assert((np.abs(vals - true_vals) < 1e-6).all()) + + def test_to_t_pdf(self): + f = self.S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100) + vals = ['%2.3f' % val for val in f.data[:10]] + truevals = ['0.000', '0.014', '0.027', '0.040', + '0.050', '0.059', '0.067', '0.073', '0.077', '0.082'] + for t, v in zip(truevals, vals): + assert(t == v) + + # estimated error bounds + vals = ['%2.4f' % val for val in f.err[:10]] + truevals = ['0.0000', '0.0003', '0.0003', '0.0004', + '0.0006', '0.0008', '0.0016', '0.0019', '0.0020', '0.0021'] + for t, v in zip(truevals, vals): + assert(t == v) + + @slow + def test_sim(self): + S = self.S + + import scipy.stats as st + x2 = S.sim(20000, 20) + truth1 = [0, np.sqrt(S.moment(1)[0]), 0., 0.] + funs = [np.mean, np.std, st.skew, st.kurtosis] + for fun, trueval in zip(funs, truth1): + res = fun(x2[:, 1::], axis=0) + m = res.mean() + sa = res.std() + assert(np.abs(m - trueval) < sa) + + @slow + def test_sim_nl(self): + S = self.S + + import scipy.stats as st + x2, _x1 = S.sim_nl(ns=20000, cases=40) + truth1 = [0, np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') + truth1[-1] = truth1[-1] - 3 + + # truth1 + # [0, 1.7495200310090633, 0.18673120577479801, 0.061988521262417606] + + funs = [np.mean, np.std, st.skew, st.kurtosis] + for fun, trueval in zip(funs, truth1): + res = fun(x2.data, axis=0) + m = res.mean() + sa = res.std() + # trueval, m, sa + assert(np.abs(m - trueval) < 2 * sa) + + def test_stats_nl(self): + S = self.S + me, va, sk, ku = S.stats_nl(moments='mvsk') + assert(me == 0.0) + assert_array_almost_equal(va, 3.0608203389019537) + assert_array_almost_equal(sk, 0.18673120577479801) + assert_array_almost_equal(ku, 3.0619885212624176) + + def test_testgaussian(self): + Hs = self.Sj.Hm0 + S0 = self.S + # ns =100; dt = .2 + # x1 = S0.sim(ns, dt=dt) + S = S0.copy() + me, _va, sk, ku = S.stats_nl(moments='mvsk') + S.tr = wtm.TrHermite( + mean=me, sigma=Hs / 4, skew=sk, kurt=ku, ysigma=Hs / 4) + ys = wo.mat2timeseries(S.sim(ns=2 ** 13)) + g0, _gemp = ys.trdata() + t0 = g0.dist2gauss() + t1 = S0.testgaussian(ns=2 ** 13, test0=None, cases=50) + assert(sum(t1 > t0) < 5) + + +class TestSpectrumHs5(unittest.TestCase): + def setUp(self): + self.Sj = sm.Jonswap(Hm0=5.0) + self.S = self.Sj.tospecdata() + + def test_moment(self): + S = self.S + vals, txt = S.moment() + true_vals = [1.5614600345079888, 0.95567089481941048] + true_txt = ['m0', 'm0tt'] + + assert_array_almost_equal(vals, true_vals) + for tv, v in zip(true_txt, txt): + assert(tv == v) + + def test_nyquist_freq(self): + S = self.S + assert_array_almost_equal(S.nyquist_freq(), 3.0) + + def test_sampling_period(self): + S = self.S + assert_array_almost_equal(S.sampling_period(), 1.0471975511965976) + + def test_normalize(self): + S = self.S + mom, txt = S.moment(2) + assert_array_almost_equal(mom, + [1.5614600345079888, 0.95567089481941048]) + assert_array_equal(txt, ['m0', 'm0tt']) + vals, _txt = S.moment(2) + true_vals = [1.5614600345079888, 0.95567089481941048] + assert_array_almost_equal(vals, true_vals) + + Sn = S.copy() + Sn.normalize() + + # Now the moments should be one + new_vals, _txt = Sn.moment(2) + assert_array_almost_equal(new_vals, np.ones(2)) + + def test_characteristic(self): + S = self.S + ch, R, txt = S.characteristic(1) + assert_array_almost_equal(ch, 8.59007646) + assert_array_almost_equal(R, 0.03040216) + self.assert_(txt == ['Tm01']) + + ch, R, txt = S.characteristic([1, 2, 3]) # fact a vector of integers + assert_array_almost_equal(ch, [8.59007646, 8.03139757, 5.62484314]) + assert_array_almost_equal(R, + [[0.03040216, 0.02834263, NAN], + [0.02834263, 0.0274645, NAN], + [NAN, NAN, 0.01500249]]) + assert_array_equal(txt, ['Tm01', 'Tm02', 'Tm24']) + + ch, R, txt = S.characteristic('Ss') # fact a string + assert_array_almost_equal(ch, [0.04963112]) + assert_array_almost_equal(R, [[2.63624782e-06]]) + assert_array_equal(txt, ['Ss']) + + # fact a list of strings + ch, R, txt = S.characteristic(['Hm0', 'Tm02']) + + assert_array_almost_equal(ch, + [4.99833578, 8.03139757]) + assert_array_almost_equal(R, [[0.05292989, 0.02511371], + [0.02511371, 0.0274645]]) + assert_array_equal(txt, ['Hm0', 'Tm02']) + + +class TestSpectrumHs3(unittest.TestCase): + def test_bandwidth(self): + + Sj = sm.Jonswap(Hm0=3, Tp=7) + w = np.linspace(0, 4, 256) + S = SpecData1D(Sj(w), w) # Make spectrum object from numerical values + vals = S.bandwidth([0, 1, 2, 3]) + true_vals = [0.73062845, 0.34476034, 0.68277527, 2.90817052] + assert_array_almost_equal(vals, true_vals) -def test_bandwidth(): - - Sj = sm.Jonswap(Hm0=3, Tp=7) - w = np.linspace(0, 4, 256) - S = SpecData1D(Sj(w), w) # Make spectrum object from numerical values - vals = S.bandwidth([0, 1, 2, 3]) - true_vals = np.array([0.73062845, 0.34476034, 0.68277527, 2.90817052]) - assert((np.abs(vals - true_vals) < 1e-7).all()) - - -def test_docstrings(): - import doctest - doctest.testmod() if __name__ == '__main__': import nose nose.run() - # test_docstrings() - # test_tocovdata() - # test_tocovmatrix() - # test_sim() - # test_bandwidth() diff --git a/wafo/stats/_distn_infrastructure.py b/wafo/stats/_distn_infrastructure.py index 0842a9d..8bc802c 100644 --- a/wafo/stats/_distn_infrastructure.py +++ b/wafo/stats/_distn_infrastructure.py @@ -1,9 +1,11 @@ -from scipy.stats._distn_infrastructure import * -from scipy.stats._distn_infrastructure import (_skew, _kurtosis, # @UnresolvedImport - _lazywhere, _ncx2_log_pdf, _ncx2_pdf, _ncx2_cdf) -from wafo.stats.estimation import FitDistribution -from wafo.stats._constants import _EPS, _XMAX -import numpy as np +from __future__ import absolute_import +from scipy.stats._distn_infrastructure import * # @UnusedWildImport +from scipy.stats._distn_infrastructure import (_skew, # @UnusedImport + _kurtosis, _lazywhere, _ncx2_log_pdf, # @IgnorePep8 @UnusedImport + _ncx2_pdf, _ncx2_cdf) # @UnusedImport @IgnorePep8 +from .estimation import FitDistribution +from ._constants import _XMAX + _doc_default_example = """\ Examples @@ -326,10 +328,11 @@ def nlogps(self, theta, x): return T -def _reduce_func(self, args, kwds): +def _reduce_func(self, args, options): # First of all, convert fshapes params to fnum: eg for stats.beta, # shapes='a, b'. To fix `a`, can specify either `f1` or `fa`. # Convert the latter into the former. + kwds = options.copy() if self.shapes: shapes = self.shapes.replace(',', ' ').split() for j, s in enumerate(shapes): @@ -384,9 +387,9 @@ def _reduce_func(self, args, kwds): return x0, func, restore, args -def fit(self, data, *args, **kwds): +def fit(self, data, *args, **kwargs): """ - Return ML or MPS estimate for shape, location, and scale parameters from data. + Return ML/MPS estimate for shape, location, and scale parameters from data. ML and MPS stands for Maximum Likelihood and Maximum Product Spacing, respectively. Starting estimates for @@ -476,6 +479,7 @@ def fit(self, data, *args, **kwds): if Narg > self.numargs: raise TypeError("Too many input arguments.") + kwds = kwargs.copy() start = [None]*2 if (Narg < self.numargs) or not ('loc' in kwds and 'scale' in kwds): @@ -573,4 +577,3 @@ rv_continuous.nlogps = nlogps rv_continuous._reduce_func = _reduce_func rv_continuous.fit = fit rv_continuous.fit2 = fit2 - diff --git a/wafo/tests/__init__.py b/wafo/tests/__init__.py index 8b13789..e69de29 100644 --- a/wafo/tests/__init__.py +++ b/wafo/tests/__init__.py @@ -1 +0,0 @@ - diff --git a/wafo/tests/conftest.py b/wafo/tests/conftest.py index 356cf9c..5c8859b 100644 --- a/wafo/tests/conftest.py +++ b/wafo/tests/conftest.py @@ -9,4 +9,4 @@ """ from __future__ import print_function, absolute_import, division -import pytest +import pytest # @UnusedImport diff --git a/wafo/tests/test_integrate.py b/wafo/tests/test_integrate.py index 4001813..f538699 100644 --- a/wafo/tests/test_integrate.py +++ b/wafo/tests/test_integrate.py @@ -7,10 +7,23 @@ import unittest import numpy as np from numpy import exp, Inf from numpy.testing import assert_array_almost_equal -from wafo.integrate import gaussq +from wafo.integrate import gaussq, quadgr, clencurt, romberg -class Gaussq(unittest.TestCase): +class TestIntegrators(unittest.TestCase): + def test_clencurt(self): + val, err = clencurt(np.exp, 0, 2) + assert_array_almost_equal(val, np.expm1(2)) + self.assert_(err < 1e-10) + + def test_romberg(self): + tol = 1e-7 + q, err = romberg(np.sqrt, 0, 10, 0, abseps=tol) + assert_array_almost_equal(q, 2.0/3 * 10**(3./2)) + self.assert_(err < tol) + + +class TestGaussq(unittest.TestCase): ''' 1 : p(x) = 1 a =-1, b = 1 Gauss-Legendre 2 : p(x) = exp(-x^2) a =-inf, b = inf Hermite @@ -60,6 +73,52 @@ class Gaussq(unittest.TestCase): val, _err = gaussq(lambda x: x, 0, 1, wfun=9) assert_array_almost_equal(val, 0.26666667) + +class TestQuadgr(unittest.TestCase): + def test_log(self): + Q, err = quadgr(np.log, 0, 1) + assert_array_almost_equal(Q, -1) + self.assert_(err < 1e-5) + + def test_exp(self): + Q, err = quadgr(np.exp, 0, 9999*1j*np.pi) + assert_array_almost_equal(Q, -2.0000000000122662) + self.assert_(err < 1.0e-8) + + def test_integral3(self): + tol = 1e-12 + Q, err = quadgr(lambda x: np.sqrt(4-x**2), 0, 2, tol) + assert_array_almost_equal(Q, np.pi) + self.assert_(err < tol) + # (3.1415926535897811, 1.5809575870662229e-13) + + def test_integral4(self): + Q, err = quadgr(lambda x: 1./x**0.75, 0, 1) + assert_array_almost_equal(Q, 4) + self.assert_(err < 1.0e-12) + + def test_integrand4(self): + tol = 1e-10 + Q, err = quadgr(lambda x: 1./np.sqrt(1-x**2), -1, 1, tol) + assert_array_almost_equal(Q, np.pi) + self.assert_(err < tol) + # (3.141596056985029, 6.2146261559092864e-06) + + def test_integrand5(self): + tol = 1e-9 + Q, err = quadgr(lambda x: np.exp(-x**2), -np.inf, np.inf, tol) + + assert_array_almost_equal(Q, np.sqrt(np.pi)) + self.assert_(err < tol) + # (1.7724538509055152, 1.9722334876348668e-11) + + def test_integrand6(self): + tol = 1e-9 + Q, err = quadgr(lambda x: np.cos(x)*np.exp(-x), 0, np.inf, tol) + assert_array_almost_equal(Q, 0.5) + self.assert_(err < tol) + # (0.50000000000000044, 7.3296813063450372e-11) + if __name__ == "__main__": # import sys;sys.argv = ['', 'Test.testName'] unittest.main() diff --git a/wafo/transform/models.py b/wafo/transform/models.py index 0d22d35..3efe47c 100644 --- a/wafo/transform/models.py +++ b/wafo/transform/models.py @@ -478,8 +478,8 @@ class TrOchi(TrCommon2): ''' Returns ga, gb, sigma2, mean2 ''' - if (self._phat is None or self.sigma != self._phat[0] - or self.mean != self._phat[1]): + if (self._phat is None or self.sigma != self._phat[0] or + self.mean != self._phat[1]): self._par_from_stats() # sigma1 = self._phat[0] # mean1 = self._phat[1] diff --git a/wafo/wave_theory/core.py b/wafo/wave_theory/core.py index e844962..297f7c5 100644 --- a/wafo/wave_theory/core.py +++ b/wafo/wave_theory/core.py @@ -7,9 +7,11 @@ from __future__ import absolute_import import numpy as np from numpy import exp, expm1, inf, nan, pi, hstack, where, atleast_1d, cos, sin from .dispersion_relation import w2k, k2w # @UnusedImport - +from ..misc import JITImport __all__ = ['w2k', 'k2w', 'sensor_typeid', 'sensor_type', 'TransferFunction'] +_MODELS = JITImport('wafo.spectrum.models') + def hyperbolic_ratio(a, b, sa, sb): ''' @@ -349,7 +351,8 @@ class TransferFunction(object): Gwt = -Gwt return Hw, Gwt __call__ = tran -#---Private member methods + +# --- Private member methods --- def _get_ee_cthxy(self, theta, kw): # convert from angle in degrees to radians @@ -358,7 +361,7 @@ class TransferFunction(object): thyr = self.thetay * pi / 180 cthx = bet * cos(theta - thxr + pi / 2) - #cthy = cos(theta-thyr-pi/2) + # cthy = cos(theta-thyr-pi/2) cthy = bet * sin(theta - thyr) # Compute location complex exponential @@ -380,14 +383,14 @@ class TransferFunction(object): zk = kw * z # z measured positive upward from sea floor return zk - #--- Surface elevation --- + # --- Surface elevation --- def _n(self, w, theta, kw): '''n = Eta = wave profile ''' ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) return np.ones_like(w), ee - #---- Vertical surface velocity and acceleration----- + # --- Vertical surface velocity and acceleration --- def _n_t(self, w, theta, kw): ''' n_t = Eta_t ''' ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -398,7 +401,7 @@ class TransferFunction(object): ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) return w ** 2, -ee - #--- Surface slopes --- + # --- Surface slopes --- def _n_x(self, w, theta, kw): ''' n_x = Eta_x = x-slope''' ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -409,7 +412,7 @@ class TransferFunction(object): ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw) return kw, 1j * cthy * ee - #--- Surface curvatures --- + # --- Surface curvatures --- def _n_xx(self, w, theta, kw): ''' n_xx = Eta_xx = Surface curvature (x-dir)''' ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -425,7 +428,7 @@ class TransferFunction(object): ee, cthx, cthy = self._get_ee_cthxy(theta, kw) return kw ** 2, -cthx * cthy * ee - #--- Pressure--- + # --- Pressure--- def _p(self, w, theta, kw): ''' pressure fluctuations''' ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -434,7 +437,7 @@ class TransferFunction(object): # hyperbolic_ratio = cosh(zk)/cosh(hk) return self.rho * self.g * hyperbolic_ratio(zk, hk, 1, 1), ee - #---- Water particle velocities --- + # --- Water particle velocities --- def _u(self, w, theta, kw): ''' U = x-velocity''' ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -459,7 +462,7 @@ class TransferFunction(object): # w*sinh(zk)/sinh(hk), -? return w * hyperbolic_ratio(zk, hk, -1, -1), -1j * ee - #---- Water particle acceleration --- + # --- Water particle acceleration --- def _u_t(self, w, theta, kw): ''' U_t = x-acceleration''' ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -484,7 +487,7 @@ class TransferFunction(object): # w*sinh(zk)/sinh(hk), ? return (w ** 2) * hyperbolic_ratio(zk, hk, -1, -1), -ee - #---- Water particle displacement --- + # --- Water particle displacement --- def _x_p(self, w, theta, kw): ''' X_p = x-displacement''' ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) @@ -508,88 +511,73 @@ class TransferFunction(object): zk = self._get_zk(kw) return hyperbolic_ratio(zk, hk, -1, -1), ee # sinh(zk)./sinh(hk), ee -# def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028): -# ''' -# Calculate pressure amplitude due to water waves. -# -# Parameters -# ---------- -# z : array-like -# depth where pressure is calculated [m] -# Hm0 : array-like -# significant wave height (same as the average of the 1/3'rd highest -# waves in a seastate. [m] -# h : real scalar -# waterdepth (default 10000 [m]) -# g : real scalar -# acceleration of gravity (default 9.81 m/s**2) -# rho : real scalar -# water density (default 1028 kg/m**3) -# -# -# Returns -# ------- -# p : ndarray -# pressure amplitude due to water waves at water depth z. [Pa] -# -# PRESSURE calculate pressure amplitude due to water waves according to -# linear theory. -# -# Example -# ----- -# >>> import pylab as plt -# >>> z = -np.linspace(10,20) -# >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20)) -# >>> plt.show() -# -# See also -# -------- -# w2k -# -# -# u = psweep.Fn*sqrt(mgf.length*9.81) -# z = -10; h = inf; -# Hm0 = 1.5;Tp = 4*sqrt(Hm0); -# S = jonswap([],[Hm0,Tp]); -# Hw = tran(S.w,0,[0 0 -z],'P',h) -# Sm = S; -# Sm.S = Hw.'.*S.S; -# x1 = spec2sdat(Sm,1000); -# pwave = pressure(z,Hm0,h) -# -# plot(psweep.x{1}/u, psweep.f) -# hold on -# plot(x1(1:100,1)-30,x1(1:100,2),'r') -# ''' -# -# -# Assume seastate with jonswap spectrum: -# -# Tp = 4 * np.sqrt(Hm0) -# gam = jonswap_peakfact(Hm0, Tp) -# Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam) -# w = 2 * np.pi / Tm02 -# kw, unused_kw2 = w2k(w, 0, h) -# -# hk = kw * h -# zk1 = kw * z -# zk = hk + zk1 # z measured positive upward from mean water level (default) -# zk = hk-zk1; % z measured positive downward from mean water level -# zk1 = -zk1; -# zk = zk1; % z measured positive upward from sea floor -# -# cosh(zk)/cosh(hk) approx exp(zk) for large h -# hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk) -# pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1)) -# pr = hyperbolic_ratio(zk, hk, 1, 1) -# pressure = (rho * g * Hm0 / 2) * pr -# -## pos = [np.zeros_like(z),np.zeros_like(z),z] -## tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g) -## Hw, Gwt = tf.tran(w,0) -## pressure2 = np.abs(Hw) * Hm0 / 2 -# -# return pressure + +def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028): + ''' + Calculate pressure amplitude due to water waves. + + Parameters + ---------- + z : array-like + depth where pressure is calculated [m] + Hm0 : array-like + significant wave height (same as the average of the 1/3'rd highest + waves in a seastate. [m] + h : real scalar + waterdepth (default 10000 [m]) + g : real scalar + acceleration of gravity (default 9.81 m/s**2) + rho : real scalar + water density (default 1028 kg/m**3) + + + Returns + ------- + p : ndarray + pressure amplitude due to water waves at water depth z. [Pa] + + PRESSURE calculate pressure amplitude due to water waves according to + linear theory. + + Example + ----- + >>> import pylab as plt + >>> z = -np.linspace(10,20) + >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20)) + >>> plt.show() + + See also + -------- + w2k + + ''' + + # Assume seastate with jonswap spectrum: + Tp = 4 * np.sqrt(Hm0) + gam = _MODELS.jonswap_peakfact(Hm0, Tp) + Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam) + w = 2 * np.pi / Tm02 + kw, unused_kw2 = w2k(w, 0, h) + + hk = kw * h + zk1 = kw * z + zk = hk + zk1 # z measured positive upward from mean water level (default) + # zk = hk-zk1 # z measured positive downward from mean water level + # zk1 = -zk1 + # zk = zk1 # z measured positive upward from sea floor + + # cosh(zk)/cosh(hk) approx exp(zk) for large h + # hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk) + # pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1)) + pr = hyperbolic_ratio(zk, hk, 1, 1) + pressure = (rho * g * Hm0 / 2) * pr + + # pos = [np.zeros_like(z),np.zeros_like(z),z] + # tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g) + # Hw, Gwt = tf.tran(w,0) + # pressure2 = np.abs(Hw) * Hm0 / 2 + + return pressure def test_docstrings():