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

Added SpecData2D.moment

Browse Source
master
per.andreas.brodtkorb 14 years ago
parent af41d6bcb6
commit a2a5a2a8be
1 changed files with 92 additions and 133 deletions
Show Stats Download Patch File Download Diff File

221
pywafo/src/wafo/spectrum/core.py
Unescape Escape View File

@ -1,5 +1,5 @@
from __future__ import division from __future__ import division
from wafo.misc import meshgrid from wafo.misc import meshgrid, gravity
from wafo.objects import mat2timeseries, TimeSeries from wafo.objects import mat2timeseries, TimeSeries
import warnings import warnings
@ -2931,7 +2931,7 @@ class SpecData2D(WafoData):
self.phi = mod(self.phi+phi0+pi,2*pi)-pi self.phi = mod(self.phi+phi0+pi,2*pi)-pi
if (rotateGrid and (self.phi!=0)): if (rotateGrid and (self.phi!=0)):
# Do a physical rotation of spectrum # Do a physical rotation of spectrum
theta = Snew.args[0] #theta = Snew.args[0]
ntOld = len(theta); ntOld = len(theta);
if (mod(theta[0]-theta[-1],2*pi)==0): if (mod(theta[0]-theta[-1],2*pi)==0):
nt = ntOld-1 nt = ntOld-1
@ -2939,7 +2939,7 @@ class SpecData2D(WafoData):
nt = ntOld nt = ntOld
theta[0:nt] = mod(theta[0:nt]-self.phi+pi,2*pi)-pi theta[0:nt] = mod(theta[0:nt]-self.phi+pi,2*pi)-pi
Snew.phi = 0; self.phi = 0
ind = theta.argsort() ind = theta.argsort()
self.data = self.data[ind,:] self.data = self.data[ind,:]
self.args[0] = theta[ind] self.args[0] = theta[ind]
@ -2980,7 +2980,7 @@ class SpecData2D(WafoData):
def moment(self, nr=2, vari='xt', even=True): def moment(self, nr=2, vari='xt'):
''' '''
Calculates spectral moments from spectrum Calculates spectral moments from spectrum
@ -2991,10 +2991,6 @@ class SpecData2D(WafoData):
vari : string vari : string
variables in model, optional when two-dim.spectrum, variables in model, optional when two-dim.spectrum,
string with 'x' and/or 'y' and/or 't' string with 'x' and/or 'y' and/or 't'
even : bool
False for all moments,
True for only even orders
Returns Returns
------- -------
m : list of moments m : list of moments
@ -3024,8 +3020,15 @@ class SpecData2D(WafoData):
For the calculation of moments see Baxevani et al. For the calculation of moments see Baxevani et al.
Example: Example:
S=demospec('dir') >>> import wafo.spectrum.models as sm
[m,mtext]=spec2mom(S,2,'xyt') >>> D = sm.Spreading()
>>> SD = D.tospecdata2d(sm.Jonswap().tospecdata())
>>> m,mtext = SD.moment(nr=2,vari='xyt')
>>> m,mtext
(array([ 3.06081401e+00, 1.31633777e-01, -1.06245690e-17,
2.13024554e+00, 1.14268485e-02, 8.07334609e-03,
1.67721091e+00, -5.29341959e-19, 1.08960751e-01,
1.08960751e-01]), ['m0', 'mx', 'my', 'mt', 'mxx', 'myy', 'mtt', 'mxy', 'mxt', 'myt'])
References References
---------- ----------
@ -3058,134 +3061,90 @@ class SpecData2D(WafoData):
w = ravel(S1.args[0]) w = ravel(S1.args[0])
theta = S1.args[1]-S1.phi theta = S1.args[1]-S1.phi
S = S1.data S = S1.data
Sw = simps(S,x=theta) Sw = simps(S,x=theta,axis=0)
m = [simps(Sw,x=w)] m = [simps(Sw,x=w)]
mtext=['m0'] mtext=['m0']
if nr>0: if nr>0:
vec = []
g = np.atleast_1d(S1.__dict__.get('g', gravity()))
kx=w**2/g[0] # maybe different normalization in x and y => diff. g
ky=w**2/g[-1]
nw=w.size nw=w.size
if 'x' in vari: if 'x' in vari:
Sc = simps(S*np.cos(theta[:,None]),x=theta) ct = np.cos(theta[:,None])
#% integral S*cos(th) dth Sc = simps(S*ct,x=theta, axis=0)
#end vec.append(kx*Sc)
mtext.append('mx')
if 'y' in vari: if 'y' in vari:
Ss = simps(S*np.sin(theta[:,None]),x=theta) st = np.sin(theta[:,None])
#end Ss = simps(S*st,x=theta, axis=0)
# if ~isfield(S1,'g') vec.append(ky*Ss)
# S1.g=gravity mtext.append('my')
# end if 't' in vari:
# kx=w.^2/S1.g(1) % maybe different normalization in x and y => diff. g vec.append(w*Sw)
# ky=w.^2/S1.g(end) mtext.append('mt')
#
# if Nv>=1 if nr>1:
# switch vari if 'x' in vari:
# case 'x' Sc2 = simps(S*ct**2,x=theta, axis=0)
# vec = kx.*Sc vec.append(kx**2*Sc2)
# mtext(end+1)={'mx'} mtext.append('mxx')
# case 'y' if 'y' in vari:
# vec = ky.*Ss Ss2 = simps(S*st**2,x=theta, axis=0)
# mtext(end+1)={'my'} vec.append(ky**2*Ss2)
# case 't' mtext.append('myy')
# vec = w.*Sw if 't' in vari:
# mtext(end+1)={'mt'} vec.append(w**2*Sw)
# end mtext.append('mtt')
# else if 'x' in vari and 'y' in vari:
# vec = [kx.*Sc ky.*Ss w*Sw] Scs = simps(S*ct*st,x=theta, axis=0)
# mtext(end+(1:3))={'mx', 'my', 'mt'} vec.append(kx*ky*Scs)
# end mtext.append('mxy')
# if nr>1 if 'x' in vari and 't' in vari:
# if strcmpi(vari(1),'x') vec.append(kx*w*Sc)
# Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).' mtext.append('mxt')
# % integral S*cos(th) dth if 'y' in vari and 't' in vari:
# Sc2=simpson(th,S1.S.*(cos(th).^2*ones(1,nw))).' vec.append(ky*w*Sc)
# % integral S*cos(th)^2 dth mtext.append('myt')
# end
# if strcmpi(vari(1),'y')||strcmpi(vari(2),'y') if nr>3:
# Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).' if 'x' in vari:
# % integral S*sin(th) dth Sc3 = simps(S*ct**3,x=theta, axis=0)
# Ss2=simpson(th,S1.S.*(sin(th).^2*ones(1,nw))).' Sc4 = simps(S*ct**4,x=theta, axis=0)
# % integral S*sin(th)^2 dth vec.append(kx**4*Sc4)
# if strcmpi(vari(1),'x') mtext.append('mxxxx')
# Scs=simpson(th,S1.S.*((cos(th).*sin(th))*ones(1,nw))).' if 'y' in vari:
# % integral S*cos(th)*sin(th) dth Ss3 = simps(S*st**3,x=theta, axis=0)
# end Ss4 = simps(S*st**4,x=theta, axis=0)
# end vec.append(ky**4*Ss4)
# if ~isfield(S1,'g') mtext.append('myyyy')
# S1.g=gravity if 't' in vari:
# end vec.append(w**4*Sw)
# mtext.append('mtttt')
# if Nv==2
# switch vari if 'x' in vari and 'y' in vari:
# case 'xy' Sc2s = simps(S*ct**2*st,x=theta, axis=0)
# vec=[kx.*Sc ky.*Ss kx.^2.*Sc2 ky.^2.*Ss2 kx.*ky.*Scs] Sc3s = simps(S*ct**3*st,x=theta, axis=0)
# mtext(end+(1:5))={'mx','my','mxx', 'myy', 'mxy'} Scs2 = simps(S*ct*st**2,x=theta, axis=0)
# case 'xt' Scs3 = simps(S*ct*st**3,x=theta, axis=0)
# vec=[kx.*Sc w.*Sw kx.^2.*Sc2 w.^2.*Sw kx.*w.*Sc] Sc2s2 = simps(S*ct**2*st**2,x=theta, axis=0)
# mtext(end+(1:5))={'mx','mt','mxx', 'mtt', 'mxt'} vec.extend((kx**3*ky*Sc3s,kx**2*ky**2*Sc2s2, kx*ky**3*Scs3))
# case 'yt' mtext.extend(('mxxxy','mxxyy','mxyyy'))
# vec=[ky.*Ss w.*Sw ky.^2.*Ss2 w.^2.*Sw ky.*w.*Ss] if 'x' in vari and 't' in vari:
# mtext(end+(1:5))={'my','mt','myy', 'mtt', 'myt'} vec.extend((kx**3*w*Sc3, kx**2*w**2*Sc2, kx*w**3*Sc))
# end mtext.extend(('mxxxt','mxxtt','mxttt'))
# else if 'y' in vari and 't' in vari:
# vec=[kx.*Sc ky.*Ss w.*Sw kx.^2.*Sc2 ky.^2.*Ss2 w.^2.*Sw kx.*ky.*Scs kx.*w.*Sc ky.*w.*Ss] vec.extend((ky**3*w*Ss3, ky**2*w**2*Ss2, ky*w**3*Ss))
# mtext(end+(1:9))={'mx','my','mt','mxx', 'myy', 'mtt', 'mxy', 'mxt', 'myt'} mtext.extend(('myyyt','myytt','myttt'))
# end if 'x' in vari and 'y' in vari and 't' in vari:
# if nr>3 vec.extend((kx**2*ky*w*Sc2s, kx*ky**2*w*Scs2, kx*ky*w**2*Scs))
# if strcmpi(vari(1),'x') mtext.extend(('mxxyt','mxyyt','mxytt'))
# Sc3=simpson(th,S1.S.*(cos(th).^3*ones(1,nw))).' #end % if nr>1
# % integral S*cos(th)^3 dth m.extend([simps(vals, x=w) for vals in vec])
# Sc4=simpson(th,S1.S.*(cos(th).^4*ones(1,nw))).' return np.asarray(m), mtext
# % integral S*cos(th)^4 dth
# end
# if strcmpi(vari(1),'y')||strcmpi(vari(2),'y')
# Ss3=simpson(th,S1.S.*(sin(th).^3*ones(1,nw))).'
# % integral S*sin(th)^3 dth
# Ss4=simpson(th,S1.S.*(sin(th).^4*ones(1,nw))).'
# % integral S*sin(th)^4 dth
# if strcmpi(vari(1),'x') %both x and y
# Sc2s=simpson(th,S1.S.*((cos(th).^2.*sin(th))*ones(1,nw))).'
# % integral S*cos(th)^2*sin(th) dth
# Sc3s=simpson(th,S1.S.*((cos(th).^3.*sin(th))*ones(1,nw))).'
# % integral S*cos(th)^3*sin(th) dth
# Scs2=simpson(th,S1.S.*((cos(th).*sin(th).^2)*ones(1,nw))).'
# % integral S*cos(th)*sin(th)^2 dth
# Scs3=simpson(th,S1.S.*((cos(th).*sin(th).^3)*ones(1,nw))).'
# % integral S*cos(th)*sin(th)^3 dth
# Sc2s2=simpson(th,S1.S.*((cos(th).^2.*sin(th).^2)*ones(1,nw))).'
# % integral S*cos(th)^2*sin(th)^2 dth
# end
# end
# if Nv==2
# switch vari
# case 'xy'
# vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 kx.^3.*ky.*Sc3s ...
# kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3]
# mtext(end+(1:5))={'mxxxx','myyyy','mxxxy','mxxyy','mxyyy'}
# case 'xt'
# vec=[vec kx.^4.*Sc4 w.^4.*Sw kx.^3.*w.*Sc3 ...
# kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc]
# mtext(end+(1:5))={'mxxxx','mtttt','mxxxt','mxxtt','mxttt'}
# case 'yt'
# vec=[vec ky.^4.*Ss4 w.^4.*Sw ky.^3.*w.*Ss3 ...
# ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss]
# mtext(end+(1:5))={'myyyy','mtttt','myyyt','myytt','myttt'}
# end
# else
# vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 w.^4.*Sw kx.^3.*ky.*Sc3s ...
# kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3 kx.^3.*w.*Sc3 ...
# kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc ky.^3.*w.*Ss3 ...
# ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss kx.^2.*ky.*w.*Sc2s ...
# kx.*ky.^2.*w.*Scs2 kx.*ky.*w.^2.*Scs]
# mtext(end+(1:15))={'mxxxx','myyyy','mtttt','mxxxy','mxxyy',...
# 'mxyyy','mxxxt','mxxtt','mxttt','myyyt','myytt','myttt','mxxyt','mxyyt','mxytt'}
#
# end % if Nv==2 ... else ...
# end % if nr>3
# end % if nr>1
# m=[m simpson(w,vec)]
# end % if nr>0
# % end %%if Nv==1... else... to be removed
# end % ... else two-dim spectrum

Write Preview
Loading…
Cancel
Save