From e4918839a981f2cf900acbd941ddcb7ce9c49ae0 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Fri, 15 Apr 2011 15:59:23 +0000 Subject: [PATCH] Translated matlab tran function into a TransferFunction class --- pywafo/src/wafo/objects.py | 459 +++++++++++++++++++++++++++- pywafo/src/wafo/stats/estimation.py | 2 +- 2 files changed, 454 insertions(+), 7 deletions(-) diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index dd58e19..f0907c2 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -2228,7 +2228,67 @@ class TimeSeries(WafoData): #wafostamp return figs -def sensortypeid(*sensortypes): +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, sign(sa), 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 @@ -2260,14 +2320,14 @@ def sensortypeid(*sensortypes): 17, 'Z_p' : Water particle displacement in z-direction from its mean position Example: - >>> sensortypeid('W','v') + >>> sensor_typeid('W','v') [11, 10] - >>> sensortypeid('rubbish') + >>> sensor_typeid('rubbish') [nan] See also -------- - sensortype + sensor_type ''' sensorid_table = dict(n=0, n_t=1, n_tt=2, n_x=3, n_y=4, n_xx=5, @@ -2280,7 +2340,7 @@ def sensortypeid(*sensortypes): -def sensortype(*sensorids): +def sensor_type(*sensorids): ''' Return sensortype name @@ -2312,7 +2372,7 @@ def sensortype(*sensorids): 17, 'Z_p' : Water particle displacement in z-direction from its mean position Example: - >>> sensortype(range(3)) + >>> sensor_type(range(3)) ('n', 'n_t', 'n_tt') See also @@ -2333,6 +2393,393 @@ def sensortype(*sensorids): #except: # raise ValueError('Input must be an integer!') +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 = 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 + +#---Private member methods + def _get_ee_cthxy(self, theta, kw): + # convert to 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 + def main0(): import wafo ts = wafo.objects.mat2timeseries(wafo.data.sea()) diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index b3f593b..f5afeb5 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -290,7 +290,7 @@ class Profile(object): phatfree = optimize.fmin(mylogfun, phatfree, args=(p_opt,) , disp=0) LLt = -mylogfun(phatfree, p_opt) if LLt>Lmax: - foundNewphat = True + #foundNewphat = True warnings.warn('The fitted parameters does not provide the optimum fit. Something wrong with fit') dL = Lmax-LLt self.alpha_cross_level -= dL