Small updates
per.andreas.brodtkorb
14 years ago
parent
d0834a14b8
commit
afd3c87e22
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@ -0,0 +1,2 @@
|
||||
from core import *
|
||||
import dispersion_relation
|
||||
@ -0,0 +1,562 @@
|
||||
'''
|
||||
Created on 3. juni 2011
|
||||
|
||||
@author: pab
|
||||
'''
|
||||
import numpy as np
|
||||
from numpy import exp, expm1, inf, nan, pi, hstack, where, atleast_1d, cos, sin
|
||||
from dispersion_relation import w2k, k2w
|
||||
|
||||
|
||||
def hyperbolic_ratio(a, b, sa, sb):
|
||||
'''
|
||||
Return ratio of hyperbolic functions
|
||||
to allow extreme variations of arguments.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
a, b : array-like
|
||||
arguments vectors of the same size
|
||||
sa, sb : scalar integers
|
||||
defining the hyperbolic function used, i.e., f(x,1)=cosh(x), f(x,-1)=sinh(x)
|
||||
|
||||
Returns
|
||||
-------
|
||||
r : ndarray
|
||||
f(a,sa)/f(b,sb), ratio of hyperbolic functions of same
|
||||
size as a and b
|
||||
Examples
|
||||
--------
|
||||
>>> x = [-2,0,2]
|
||||
>>> hyperbolic_ratio(x,1,1,1) # gives r=cosh(x)/cosh(1)
|
||||
array([ 2.438107 , 0.64805427, 2.438107 ])
|
||||
>>> hyperbolic_ratio(x,1,1,-1) # gives r=cosh(x)/sinh(1)
|
||||
array([ 3.20132052, 0.85091813, 3.20132052])
|
||||
>>> hyperbolic_ratio(x,1,-1,1) # gives r=sinh(x)/cosh(1)
|
||||
array([-2.35040239, 0. , 2.35040239])
|
||||
>>> hyperbolic_ratio(x,1,-1,-1) # gives r=sinh(x)/sinh(1)
|
||||
array([-3.08616127, 0. , 3.08616127])
|
||||
>>> hyperbolic_ratio(1,x,1,1) # gives r=cosh(1)/cosh(x)
|
||||
array([ 0.41015427, 1.54308063, 0.41015427])
|
||||
>>> hyperbolic_ratio(1,x,1,-1) # gives r=cosh(1)/sinh(x)
|
||||
array([-0.42545906, inf, 0.42545906])
|
||||
>>> hyperbolic_ratio(1,x,-1,1) # gives r=sinh(1)/cosh(x)
|
||||
array([ 0.3123711 , 1.17520119, 0.3123711 ])
|
||||
>>> hyperbolic_ratio(1,x,-1,-1) # gives r=sinh(1)/sinh(x)
|
||||
array([-0.32402714, inf, 0.32402714])
|
||||
|
||||
See also
|
||||
--------
|
||||
tran
|
||||
'''
|
||||
|
||||
ak, bk, sak, sbk = np.atleast_1d(a, b, np.sign(sa), np.sign(sb))
|
||||
# old call
|
||||
#return exp(ak-bk)*(1+sak*exp(-2*ak))/(1+sbk*exp(-2*bk))
|
||||
# TODO: Does not always handle division by zero correctly
|
||||
|
||||
signRatio = np.where(sak * ak < 0, sak, 1)
|
||||
signRatio = np.where(sbk * bk < 0, sbk * signRatio, signRatio)
|
||||
|
||||
bk = np.abs(bk)
|
||||
ak = np.abs(ak)
|
||||
|
||||
num = np.where(sak < 0, expm1(-2 * ak), 1 + exp(-2 * ak))
|
||||
den = np.where(sbk < 0, expm1(-2 * bk), 1 + exp(-2 * bk))
|
||||
iden = np.ones(den.shape) * inf
|
||||
ind = np.flatnonzero(den != 0)
|
||||
iden.flat[ind] = 1.0 / den[ind]
|
||||
val = np.where(num == den, 1, num * iden)
|
||||
return signRatio * exp(ak - bk) * val #((sak+exp(-2*ak))/(sbk+exp(-2*bk)))
|
||||
|
||||
def sensor_typeid(*sensortypes):
|
||||
''' Return ID for sensortype name
|
||||
|
||||
Parameter
|
||||
---------
|
||||
sensortypes : list of strings defining the sensortype
|
||||
|
||||
Returns
|
||||
-------
|
||||
sensorids : list of integers defining the sensortype
|
||||
|
||||
Valid senor-ids and -types for time series are as follows:
|
||||
0, 'n' : Surface elevation (n=Eta)
|
||||
1, 'n_t' : Vertical surface velocity
|
||||
2, 'n_tt' : Vertical surface acceleration
|
||||
3, 'n_x' : Surface slope in x-direction
|
||||
4, 'n_y' : Surface slope in y-direction
|
||||
5, 'n_xx' : Surface curvature in x-direction
|
||||
6, 'n_yy' : Surface curvature in y-direction
|
||||
7, 'n_xy' : Surface curvature in xy-direction
|
||||
8, 'P' : Pressure fluctuation about static MWL pressure
|
||||
9, 'U' : Water particle velocity in x-direction
|
||||
10, 'V' : Water particle velocity in y-direction
|
||||
11, 'W' : Water particle velocity in z-direction
|
||||
12, 'U_t' : Water particle acceleration in x-direction
|
||||
13, 'V_t' : Water particle acceleration in y-direction
|
||||
14, 'W_t' : Water particle acceleration in z-direction
|
||||
15, 'X_p' : Water particle displacement in x-direction from its mean position
|
||||
16, 'Y_p' : Water particle displacement in y-direction from its mean position
|
||||
17, 'Z_p' : Water particle displacement in z-direction from its mean position
|
||||
|
||||
Example:
|
||||
>>> sensor_typeid('W','v')
|
||||
[11, 10]
|
||||
>>> sensor_typeid('rubbish')
|
||||
[nan]
|
||||
|
||||
See also
|
||||
--------
|
||||
sensor_type
|
||||
'''
|
||||
|
||||
sensorid_table = dict(n=0, n_t=1, n_tt=2, n_x=3, n_y=4, n_xx=5,
|
||||
n_yy=6, n_xy=7, p=8, u=9, v=10, w=11, u_t=12,
|
||||
v_t=13, w_t=14, x_p=15, y_p=16, z_p=17)
|
||||
try:
|
||||
return [sensorid_table.get(name.lower(), nan) for name in sensortypes]
|
||||
except:
|
||||
raise ValueError('Input must be a string!')
|
||||
|
||||
|
||||
|
||||
def sensor_type(*sensorids):
|
||||
'''
|
||||
Return sensortype name
|
||||
|
||||
Parameter
|
||||
---------
|
||||
sensorids : vector or list of integers defining the sensortype
|
||||
|
||||
Returns
|
||||
-------
|
||||
sensornames : tuple of strings defining the sensortype
|
||||
Valid senor-ids and -types for time series are as follows:
|
||||
0, 'n' : Surface elevation (n=Eta)
|
||||
1, 'n_t' : Vertical surface velocity
|
||||
2, 'n_tt' : Vertical surface acceleration
|
||||
3, 'n_x' : Surface slope in x-direction
|
||||
4, 'n_y' : Surface slope in y-direction
|
||||
5, 'n_xx' : Surface curvature in x-direction
|
||||
6, 'n_yy' : Surface curvature in y-direction
|
||||
7, 'n_xy' : Surface curvature in xy-direction
|
||||
8, 'P' : Pressure fluctuation about static MWL pressure
|
||||
9, 'U' : Water particle velocity in x-direction
|
||||
10, 'V' : Water particle velocity in y-direction
|
||||
11, 'W' : Water particle velocity in z-direction
|
||||
12, 'U_t' : Water particle acceleration in x-direction
|
||||
13, 'V_t' : Water particle acceleration in y-direction
|
||||
14, 'W_t' : Water particle acceleration in z-direction
|
||||
15, 'X_p' : Water particle displacement in x-direction from its mean position
|
||||
16, 'Y_p' : Water particle displacement in y-direction from its mean position
|
||||
17, 'Z_p' : Water particle displacement in z-direction from its mean position
|
||||
|
||||
Example:
|
||||
>>> sensor_type(range(3))
|
||||
('n', 'n_t', 'n_tt')
|
||||
|
||||
See also
|
||||
--------
|
||||
sensor_typeid, tran
|
||||
'''
|
||||
valid_names = ('n', 'n_t', 'n_tt', 'n_x', 'n_y', 'n_xx', 'n_yy', 'n_xy',
|
||||
'p', 'u', 'v', 'w', 'u_t', 'v_t', 'w_t', 'x_p', 'y_p', 'z_p',
|
||||
nan)
|
||||
ids = atleast_1d(*sensorids)
|
||||
if isinstance(ids, list):
|
||||
ids = hstack(ids)
|
||||
n = len(valid_names) - 1
|
||||
ids = where(((ids < 0) | (n < ids)), n , ids)
|
||||
return tuple(valid_names[i] for i in ids)
|
||||
|
||||
class TransferFunction(object):
|
||||
'''
|
||||
Class for computing transfer functions based on linear wave theory
|
||||
of the system with input surface elevation,
|
||||
eta(x0,y0,t) = exp(i*(kx*x0+ky*y0-w*t)),
|
||||
and output Y determined by sensortype and position of sensor.
|
||||
|
||||
Member methods
|
||||
--------------
|
||||
tran(w, theta, kw)
|
||||
|
||||
Hw = a function of frequency only (not direction) size 1 x Nf
|
||||
Gwt = a function of frequency and direction size Nt x Nf
|
||||
w = vector of angular frequencies in Rad/sec. Length Nf
|
||||
theta = vector of directions in radians Length Nt (default 0)
|
||||
( theta = 0 -> positive x axis theta = pi/2 -> positive y axis)
|
||||
Member variables
|
||||
----------------
|
||||
pos : [x,y,z]
|
||||
vector giving coordinate position relative to [x0 y0 z0] (default [0,0,0])
|
||||
sensortype = string
|
||||
defining the sensortype or transfer function in output.
|
||||
0, 'n' : Surface elevation (n=Eta) (default)
|
||||
1, 'n_t' : Vertical surface velocity
|
||||
2, 'n_tt' : Vertical surface acceleration
|
||||
3, 'n_x' : Surface slope in x-direction
|
||||
4, 'n_y' : Surface slope in y-direction
|
||||
5, 'n_xx' : Surface curvature in x-direction
|
||||
6, 'n_yy' : Surface curvature in y-direction
|
||||
7, 'n_xy' : Surface curvature in xy-direction
|
||||
8, 'P' : Pressure fluctuation about static MWL pressure
|
||||
9, 'U' : Water particle velocity in x-direction
|
||||
10, 'V' : Water particle velocity in y-direction
|
||||
11, 'W' : Water particle velocity in z-direction
|
||||
12, 'U_t' : Water particle acceleration in x-direction
|
||||
13, 'V_t' : Water particle acceleration in y-direction
|
||||
14, 'W_t' : Water particle acceleration in z-direction
|
||||
15, 'X_p' : Water particle displacement in x-direction from its mean position
|
||||
16, 'Y_p' : Water particle displacement in y-direction from its mean position
|
||||
17, 'Z_p' : Water particle displacement in z-direction from its mean position
|
||||
h : real scalar
|
||||
water depth (default inf)
|
||||
g : real scalar
|
||||
acceleration of gravity (default 9.81 m/s**2)
|
||||
rho : real scalar
|
||||
water density (default 1028 kg/m**3)
|
||||
bet : 1 or -1
|
||||
1, theta given in terms of directions toward which waves travel (default)
|
||||
-1, theta given in terms of directions from which waves come
|
||||
igam : 1,2 or 3
|
||||
1, if z is measured positive upward from mean water level (default)
|
||||
2, if z is measured positive downward from mean water level
|
||||
3, if z is measured positive upward from sea floor
|
||||
thetax, thetay : real scalars
|
||||
angle in degrees clockwise from true north to positive x-axis and
|
||||
positive y-axis, respectively. (default theatx=90, thetay=0)
|
||||
|
||||
Example
|
||||
-------
|
||||
>>> import pylab as plt
|
||||
>>> N=50; f0=0.1; th0=0; h=50; w0 = 2*pi*f0
|
||||
>>> t = np.linspace(0,15,N)
|
||||
>>> eta0 = np.exp(-1j*w0*t)
|
||||
>>> stypes = ['n', 'n_x', 'n_y'];
|
||||
>>> tf = TransferFunction(pos=(0, 0, 0), h=50)
|
||||
>>> vals = []
|
||||
>>> fh = plt.plot(t, eta0.real, 'r.')
|
||||
>>> plt.hold(True)
|
||||
>>> for i,stype in enumerate(stypes):
|
||||
... tf.sensortype = stype
|
||||
... Hw, Gwt = tf.tran(w0,th0)
|
||||
... vals.append((Hw*Gwt*eta0).real.ravel())
|
||||
... vals[i]
|
||||
... fh = plt.plot(t, vals[i])
|
||||
>>> plt.show()
|
||||
|
||||
|
||||
See also
|
||||
--------
|
||||
dat2dspec, sensor_type, sensor_typeid
|
||||
|
||||
Reference
|
||||
---------
|
||||
Young I.R. (1994)
|
||||
"On the measurement of directional spectra",
|
||||
Applied Ocean Research, Vol 16, pp 283-294
|
||||
'''
|
||||
def __init__(self, pos=(0, 0, 0), sensortype='n', h=inf, g=9.81, rho=1028,
|
||||
bet=1, igam=1, thetax=90, thetay=0):
|
||||
self.pos = pos
|
||||
self.sensortype = sensortype if isinstance(sensortype, str) else sensor_type(sensortype)
|
||||
self.h = h
|
||||
self.g = g
|
||||
self.rho = rho
|
||||
self.bet = bet
|
||||
self.igam = igam
|
||||
self.thetax = thetax
|
||||
self.thetay = thetay
|
||||
self._tran_dict = dict(n=self._n, n_t=self._n_t, n_tt=self._n_tt,
|
||||
n_x=self._n_x, n_y=self._n_y, n_xx=self._n_xx,
|
||||
n_yy=self._n_yy, n_xy=self._n_xy,
|
||||
P=self._p, p=self._p,
|
||||
U=self._u, u=self._u,
|
||||
V=self._v, v=self._v,
|
||||
W=self._w, w=self._w,
|
||||
U_t=self._u_t, u_t=self._u_t,
|
||||
V_t=self._v_t, v_t=self._v_t,
|
||||
W_t=self._w_t, w_t=self._w_t,
|
||||
X_p=self._x_p, x_p=self._x_p,
|
||||
Y_p=self._y_p, y_p=self._y_p,
|
||||
Z_p=self._z_p, z_p=self._z_p)
|
||||
|
||||
def tran(self, w, theta=0, kw=None):
|
||||
'''
|
||||
Return transfer functions based on linear wave theory
|
||||
of the system with input surface elevation,
|
||||
eta(x0,y0,t) = exp(i*(kx*x0+ky*y0-w*t)),
|
||||
and output,
|
||||
Y = Hw*Gwt*eta, determined by sensortype and position of sensor.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
w : array-like
|
||||
vector of angular frequencies in Rad/sec. Length Nf
|
||||
theta : array-like
|
||||
vector of directions in radians Length Nt (default 0)
|
||||
( theta = 0 -> positive x axis theta = pi/2 -> positive y axis)
|
||||
kw : array-like
|
||||
vector of wave numbers corresponding to angular frequencies, w. Length Nf
|
||||
(default calculated with w2k)
|
||||
|
||||
Returns
|
||||
-------
|
||||
Hw = transfer function of frequency only (not direction) size 1 x Nf
|
||||
Gwt = transfer function of frequency and direction size Nt x Nf
|
||||
|
||||
The complete transfer function Hwt = Hw*Gwt is a function of
|
||||
w (columns) and theta (rows) size Nt x Nf
|
||||
'''
|
||||
if kw is None:
|
||||
kw, unusedkw2 = w2k(w, 0, self.h) #wave number as function of angular frequency
|
||||
|
||||
w, theta, kw = np.atleast_1d(w, theta, kw)
|
||||
# make sure they have the correct orientation
|
||||
theta.shape = (-1, 1)
|
||||
kw.shape = (-1,)
|
||||
w.shape = (-1,)
|
||||
|
||||
tran_fun = self._tran_dict[self.sensortype]
|
||||
Hw, Gwt = tran_fun(w, theta, kw)
|
||||
|
||||
# New call to avoid singularities. pab 07.11.2000
|
||||
# Set Hw to 0 for expressions w*hyperbolic_ratio(z*k,h*k,1,-1)= 0*inf
|
||||
ind = np.flatnonzero(1 - np.isfinite(Hw))
|
||||
Hw.flat[ind] = 0
|
||||
|
||||
sgn = np.sign(Hw);
|
||||
k0 = np.flatnonzero(sgn < 0)
|
||||
if len(k0): # make sure Hw>=0 ie. transfer negative signs to Gwt
|
||||
Gwt[:, k0] = -Gwt[:, k0]
|
||||
Hw[:, k0] = -Hw[:, k0]
|
||||
|
||||
if self.igam == 2:
|
||||
#pab 09 Oct.2002: bug fix
|
||||
# Changing igam by 2 should affect the directional result in the same way that changing eta by -eta!
|
||||
Gwt = -Gwt
|
||||
return Hw, Gwt
|
||||
__call__ = tran
|
||||
#---Private member methods
|
||||
def _get_ee_cthxy(self, theta, kw):
|
||||
# convert from angle in degrees to radians
|
||||
bet = self.bet
|
||||
thxr = self.thetax * pi / 180
|
||||
thyr = self.thetay * pi / 180
|
||||
|
||||
cthx = bet * cos(theta - thxr + pi / 2)
|
||||
#cthy = cos(theta-thyr-pi/2)
|
||||
cthy = bet * sin(theta - thyr)
|
||||
|
||||
# Compute location complex exponential
|
||||
x, y, unused_z = list(self.pos)
|
||||
ee = exp((1j * (x * cthx + y * cthy)) * kw) # exp(i*k(w)*(x*cos(theta)+y*sin(theta)) size Nt X Nf
|
||||
return ee, cthx, cthy
|
||||
|
||||
def _get_zk(self, kw):
|
||||
h = self.h
|
||||
z = self.pos[2]
|
||||
if self.igam == 1:
|
||||
zk = kw * (h + z) # z measured positive upward from mean water level (default)
|
||||
elif self.igam == 2:
|
||||
zk = kw * (h - z) # z measured positive downward from mean water level
|
||||
else:
|
||||
zk = kw * z # z measured positive upward from sea floor
|
||||
return zk
|
||||
|
||||
#--- 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-----
|
||||
def _n_t(self, w, theta, kw):
|
||||
''' n_t = Eta_t '''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
return w, -1j * ee;
|
||||
def _n_tt(self, w, theta, kw):
|
||||
'''n_tt = Eta_tt'''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
return w ** 2, -ee
|
||||
|
||||
#--- 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)
|
||||
return kw, 1j * cthx * ee
|
||||
def _n_y(self, w, theta, kw):
|
||||
''' n_y = Eta_y = y-slope'''
|
||||
ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
return kw, 1j * cthy * ee
|
||||
|
||||
#--- 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)
|
||||
return kw ** 2, -(cthx ** 2) * ee
|
||||
def _n_yy(self, w, theta, kw):
|
||||
''' n_yy = Eta_yy = Surface curvature (y-dir)'''
|
||||
ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
return kw ** 2, -cthy ** 2 * ee
|
||||
def _n_xy(self, w, theta, kw):
|
||||
''' n_xy = Eta_xy = Surface curvature (xy-dir)'''
|
||||
ee, cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
return kw ** 2, -cthx * cthy * ee
|
||||
|
||||
#--- Pressure---
|
||||
def _p(self, w, theta, kw):
|
||||
''' pressure fluctuations'''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return self.rho * self.g * hyperbolic_ratio(zk, hk, 1, 1), ee #hyperbolic_ratio = cosh(zk)/cosh(hk)
|
||||
|
||||
#---- Water particle velocities ---
|
||||
def _u(self, w, theta, kw):
|
||||
''' U = x-velocity'''
|
||||
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return w * hyperbolic_ratio(zk, hk, 1, -1), cthx * ee# w*cosh(zk)/sinh(hk), cos(theta)*ee
|
||||
def _v(self, w, theta, kw):
|
||||
'''V = y-velocity'''
|
||||
ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return w * hyperbolic_ratio(zk, hk, 1, -1), cthy * ee # w*cosh(zk)/sinh(hk), sin(theta)*ee
|
||||
def _w(self, w, theta, kw):
|
||||
''' W = z-velocity'''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return w * hyperbolic_ratio(zk, hk, -1, -1), -1j * ee # w*sinh(zk)/sinh(hk), -?
|
||||
|
||||
#---- Water particle acceleration ---
|
||||
def _u_t(self, w, theta, kw):
|
||||
''' U_t = x-acceleration'''
|
||||
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return (w ** 2) * hyperbolic_ratio(zk, hk, 1, -1), -1j * cthx * ee # w^2*cosh(zk)/sinh(hk), ?
|
||||
|
||||
def _v_t(self, w, theta, kw):
|
||||
''' V_t = y-acceleration'''
|
||||
ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return (w ** 2) * hyperbolic_ratio(zk, hk, 1, -1), -1j * cthy * ee # w^2*cosh(zk)/sinh(hk), ?
|
||||
def _w_t(self, w, theta, kw):
|
||||
''' W_t = z-acceleration'''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return (w ** 2) * hyperbolic_ratio(zk, hk, -1, -1), -ee # w*sinh(zk)/sinh(hk), ?
|
||||
|
||||
#---- Water particle displacement ---
|
||||
def _x_p(self, w, theta, kw):
|
||||
''' X_p = x-displacement'''
|
||||
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return hyperbolic_ratio(zk, hk, 1, -1), 1j * cthx * ee # cosh(zk)./sinh(hk), ?
|
||||
def _y_p(self, w, theta, kw):
|
||||
''' Y_p = y-displacement'''
|
||||
ee, unused_cthx, cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
zk = self._get_zk(kw)
|
||||
return hyperbolic_ratio(zk, hk, 1, -1), 1j * cthy * ee # cosh(zk)./sinh(hk), ?
|
||||
def _z_p(self, w, theta, kw):
|
||||
''' Z_p = z-displacement'''
|
||||
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
|
||||
hk = kw * self.h
|
||||
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
|
||||
|
||||
if __name__ == '__main__':
|
||||
pass
|
||||
@ -0,0 +1,205 @@
|
||||
"""
|
||||
Dispersion relation module
|
||||
--------------------------
|
||||
k2w - Translates from wave number to frequency
|
||||
w2k - Translates from frequency to wave number
|
||||
"""
|
||||
import warnings
|
||||
#import numpy as np
|
||||
from numpy import (atleast_1d, sqrt, zeros_like, arctan2, where, tanh, any, #@UnresolvedImport
|
||||
sin, cos, sign, inf, flatnonzero, finfo, cosh, abs) #@UnresolvedImport
|
||||
|
||||
__all__ = ['k2w', 'w2k']
|
||||
|
||||
def k2w(k1, k2=0e0, h=inf, g=9.81, u1=0e0, u2=0e0):
|
||||
''' Translates from wave number to frequency
|
||||
using the dispersion relation
|
||||
|
||||
Parameters
|
||||
----------
|
||||
k1 : array-like
|
||||
wave numbers [rad/m].
|
||||
k2 : array-like, optional
|
||||
second dimension wave number
|
||||
h : real scalar, optional
|
||||
water depth [m].
|
||||
g : real scalar, optional
|
||||
acceleration of gravity, see gravity
|
||||
u1, u2 : real scalars, optional
|
||||
current velocity [m/s] along dimension 1 and 2.
|
||||
note: when u1!=0 | u2!=0 then theta is not calculated correctly
|
||||
|
||||
Returns
|
||||
-------
|
||||
w : ndarray
|
||||
angular frequency [rad/s].
|
||||
theta : ndarray
|
||||
direction [rad].
|
||||
|
||||
Dispersion relation
|
||||
-------------------
|
||||
w = sqrt(g*K*tanh(K*h)) ( 0 < w < inf)
|
||||
theta = arctan2(k2,k1) (-pi < theta < pi)
|
||||
where
|
||||
K = sqrt(k1**2+k2**2)
|
||||
|
||||
The shape of w and theta is the common shape of k1 and k2 according to the
|
||||
numpy broadcasting rules.
|
||||
|
||||
See also
|
||||
--------
|
||||
w2k
|
||||
|
||||
Example
|
||||
-------
|
||||
>>> from numpy import arange
|
||||
>>> import wafo.spectrum.dispersion_relation as wsd
|
||||
>>> wsd.k2w(arange(0.01,.5,0.2))[0]
|
||||
array([ 0.3132092 , 1.43530485, 2.00551739])
|
||||
>>> wsd.k2w(arange(0.01,.5,0.2),h=20)[0]
|
||||
array([ 0.13914927, 1.43498213, 2.00551724])
|
||||
'''
|
||||
|
||||
k1i, k2i, hi, gi, u1i, u2i = atleast_1d(k1, k2, h, g, u1, u2)
|
||||
|
||||
if k1i.size == 0:
|
||||
return zeros_like(k1i)
|
||||
ku1 = k1i*u1i
|
||||
ku2 = k2i*u2i
|
||||
|
||||
theta = arctan2(k2, k1)
|
||||
|
||||
k = sqrt(k1i**2+k2i**2)
|
||||
w = where(k>0, ku1+ku2+sqrt(gi*k*tanh(k*hi)), 0.0)
|
||||
|
||||
cond = (w<0)
|
||||
if any(cond):
|
||||
txt0 = '''
|
||||
Waves and current are in opposite directions
|
||||
making some of the frequencies negative.
|
||||
Here we are forcing the negative frequencies to zero.
|
||||
'''
|
||||
warnings.warn(txt0)
|
||||
w = where(cond, 0.0, w) # force w to zero
|
||||
|
||||
return w, theta
|
||||
|
||||
def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100):
|
||||
'''
|
||||
Translates from frequency to wave number
|
||||
using the dispersion relation
|
||||
|
||||
Parameters
|
||||
----------
|
||||
w : array-like
|
||||
angular frequency [rad/s].
|
||||
theta : array-like, optional
|
||||
direction [rad].
|
||||
h : real scalar, optional
|
||||
water depth [m].
|
||||
g : real scalar or array-like of size 2.
|
||||
constant of gravity [m/s**2] or 3D normalizing constant
|
||||
|
||||
Returns
|
||||
-------
|
||||
k1, k2 : ndarray
|
||||
wave numbers [rad/m] along dimension 1 and 2.
|
||||
|
||||
Description
|
||||
-----------
|
||||
Uses Newton Raphson method to find the wave number k in the dispersion relation
|
||||
w**2= g*k*tanh(k*h).
|
||||
The solution k(w) => k1 = k(w)*cos(theta)
|
||||
k2 = k(w)*sin(theta)
|
||||
The size of k1,k2 is the common shape of w and theta according to numpy
|
||||
broadcasting rules. If w or theta is scalar it functions as a constant
|
||||
matrix of the same shape as the other.
|
||||
|
||||
Example
|
||||
-------
|
||||
>>> import pylab as plb
|
||||
>>> import wafo.spectrum.dispersion_relation as wsd
|
||||
>>> w = plb.linspace(0,3);
|
||||
>>> h = plb.plot(w,w2k(w)[0])
|
||||
>>> wsd.w2k(range(4))[0]
|
||||
array([ 0. , 0.1019368 , 0.4077472 , 0.91743119])
|
||||
>>> wsd.w2k(range(4),h=20)[0]
|
||||
array([ 0. , 0.10503601, 0.40774726, 0.91743119])
|
||||
|
||||
>>> plb.close('all')
|
||||
|
||||
See also
|
||||
--------
|
||||
k2w
|
||||
'''
|
||||
wi, th, gi = atleast_1d(w, theta, g)
|
||||
|
||||
if wi.size == 0:
|
||||
return zeros_like(wi)
|
||||
|
||||
k = 1.0*sign(wi)*wi**2.0 #% deep water
|
||||
if h > 10. ** 25:
|
||||
k2 = k*sin(th)/gi[-1] #%size np x nf
|
||||
k1 = k*cos(th)/gi[0]
|
||||
return k1, k2
|
||||
|
||||
|
||||
if gi.size > 1:
|
||||
txt0 = '''
|
||||
Finite depth in combination with 3D normalization (len(g)=2) is not implemented yet.
|
||||
'''
|
||||
raise ValueError(txt0)
|
||||
|
||||
|
||||
find = flatnonzero
|
||||
eps = finfo(float).eps
|
||||
|
||||
oshape = k.shape
|
||||
wi, k = wi.ravel(), k.ravel()
|
||||
|
||||
# Newton's Method
|
||||
# Permit no more than count_limit iterations.
|
||||
|
||||
hn = zeros_like(k)
|
||||
ix = find((wi<0) | (0<wi))
|
||||
|
||||
# Break out of the iteration loop for three reasons:
|
||||
# 1) the last update is very small (compared to x)
|
||||
# 2) the last update is very small (compared to sqrt(eps))
|
||||
# 3) There are more than 100 iterations. This should NEVER happen.
|
||||
count = 0
|
||||
while (ix.size>0 and count < count_limit):
|
||||
ki = k[ix]
|
||||
hn[ix] = (ki*tanh(ki*h)-wi[ix]**2.0/gi)/(tanh(ki*h)+ki*h/(cosh(ki*h)**2.0))
|
||||
knew = ki - hn[ix]
|
||||
# Make sure that the current guess is not zero.
|
||||
# When Newton's Method suggests steps that lead to zero guesses
|
||||
# take a step 9/10ths of the way to zero:
|
||||
ksmall = find(abs(knew)==0)
|
||||
if ksmall.size>0:
|
||||
knew[ksmall] = ki[ksmall] / 10.0
|
||||
hn[ix[ksmall]] = ki[ksmall]-knew[ksmall]
|
||||
|
||||
k[ix] = knew
|
||||
# disp(['Iteration ',num2str(count),' Number of points left: ' num2str(length(ix)) ]),
|
||||
|
||||
ix = find((abs(hn) > sqrt(eps)*abs(k)) * abs(hn) > sqrt(eps))
|
||||
count += 1
|
||||
|
||||
if count == count_limit:
|
||||
txt1 = ''' W2K did not converge.
|
||||
The maximum error in the last step was: %13.8f''' % max(hn[ix])
|
||||
warnings.warn(txt1)
|
||||
|
||||
k.shape = oshape
|
||||
|
||||
k2 = k*sin(th)
|
||||
k1 = k*cos(th)
|
||||
return k1, k2
|
||||
|
||||
def main():
|
||||
import doctest
|
||||
doctest.testmod()
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
Loading…
Reference in New Issue