TRAN Computes 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 type and pos. CALL: [Hw Gwt kw Hwt] = tran(w,theta,pos,type,h,g,rho,bet,igam,thx,thy,kw); Hwt = Hw(ones(Nt,1),:).*Gwt matrix of transfer functions values as function of w (columns) and theta (rows) size Nt x Nf 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) pos = [x,y,z] = vector giving coordinate position relative to [x0 y0 z0] (default [0,0,0]) type = Number or string defining the transfer function in output. 1, 'n' : Surface elevation (n=Eta) (default) 2, 'n_t' : Vertical surface velocity 3, 'n_tt' : Vertical surface acceleration 4, 'n_x' : Surface slope in x-direction 5, 'n_y' : Surface slope in y-direction 6, 'n_xx' : Surface curvature in x-direction 7, 'n_yy' : Surface curvature in y-direction 8, 'n_xy' : Surface curvature in xy-direction 9, 'P' : Pressure fluctuation about static MWL pressure 10, 'U' : Water particle velocity in x-direction 11, 'V' : Water particle velocity in y-direction 12, 'W' : Water particle velocity in z-direction 13, 'U_t' : Water particle acceleration in x-direction 14, 'V_t' : Water particle acceleration in y-direction 15, 'W_t' : Water particle acceleration in z-direction 16, 'X_p' : Water particle displacement in x-direction from its mean position 17, 'Y_p' : Water particle displacement in y-direction from its mean position 18, 'Z_p' : Water particle displacement in z-direction from its mean position 19, 'Dummy': Transfer function is zero h = water depth (default inf) g = acceleration of gravity (default see gravity) rho = water density (default see wdensity) bet = 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, 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 thx = angle clockwise from true north to positive x-axis in degrees (default 90) thy = angle clockwise from true north to positive y-axis in degrees (default 0) kw = vector of wave numbers corresponding to angular frequencies, w (default calculated with w2k) Example: N=5000;dt=0.4;f0=0.1;th0=0;h=50; w0 = 2*pi*f0; t = linspace(0,1500,N)'; types = [1 4 5]; bfs = [1 0 0]; pos = zeros(3,3); eta0 = exp(-i*w0*t); [Hw Gwt] = tran(w0,th0,pos(1,:),'n',h); eta = real(Hw*Gwt*eta0)+0.001*rand(N,1); [Hw Gwt] = tran(w0,th0,pos(2,:),'n_x',h); n_x = real(Hw*Gwt*eta0)+0.001*rand(N,1); [Hw Gwt] = tran(w0,th0,pos(3,:),'n_y',h); n_y = real(Hw*Gwt*eta0)+0.001*rand(N,1); S = dat2dspec([t eta n_x n_y],[pos types',bfs'],h,256,101); wspecplot(S) See also dat2dspec, wdensity, gravity
returns the constant acceleration of gravity | |
compute ratio of hyperbolic functions | |
Return ID for sensortype name | |
Translates from frequency to wave number | |
Returns the water density | |
Display message and abort function. | |
True for character array (string). | |
Write formatted data to string. |
Estimates the directional wave spectrum from timeseries | |
Creates a test case for measurement time series |
001 function [Hw, Gwt, kw,Hwt,ee]=tran(w,theta,pos,def,h,g,rho,bet,igam,thx,thy,kw); 002 %TRAN Computes transfer functions based on linear wave theory 003 % of the system with input surface elevation, 004 % eta(x0,y0,t) = exp(i*(kx*x0+ky*y0-w*t)), 005 % and output Y determined by type and pos. 006 % 007 % CALL: [Hw Gwt kw Hwt] = tran(w,theta,pos,type,h,g,rho,bet,igam,thx,thy,kw); 008 % 009 % Hwt = Hw(ones(Nt,1),:).*Gwt matrix of transfer functions values as function of 010 % w (columns) and theta (rows) size Nt x Nf 011 % Hw = a function of frequency only (not direction) size 1 x Nf 012 % Gwt = a function of frequency and direction size Nt x Nf 013 % w = vector of angular frequencies in Rad/sec. Length Nf 014 % theta = vector of directions in radians Length Nt (default 0) 015 % ( theta = 0 -> positive x axis theta = pi/2 -> positive y axis) 016 % pos = [x,y,z] = vector giving coordinate position relative to [x0 y0 z0] (default [0,0,0]) 017 % type = Number or string defining the transfer function in output. 018 % 1, 'n' : Surface elevation (n=Eta) (default) 019 % 2, 'n_t' : Vertical surface velocity 020 % 3, 'n_tt' : Vertical surface acceleration 021 % 4, 'n_x' : Surface slope in x-direction 022 % 5, 'n_y' : Surface slope in y-direction 023 % 6, 'n_xx' : Surface curvature in x-direction 024 % 7, 'n_yy' : Surface curvature in y-direction 025 % 8, 'n_xy' : Surface curvature in xy-direction 026 % 9, 'P' : Pressure fluctuation about static MWL pressure 027 % 10, 'U' : Water particle velocity in x-direction 028 % 11, 'V' : Water particle velocity in y-direction 029 % 12, 'W' : Water particle velocity in z-direction 030 % 13, 'U_t' : Water particle acceleration in x-direction 031 % 14, 'V_t' : Water particle acceleration in y-direction 032 % 15, 'W_t' : Water particle acceleration in z-direction 033 % 16, 'X_p' : Water particle displacement in x-direction from its mean position 034 % 17, 'Y_p' : Water particle displacement in y-direction from its mean position 035 % 18, 'Z_p' : Water particle displacement in z-direction from its mean position 036 % 19, 'Dummy': Transfer function is zero 037 % h = water depth (default inf) 038 % g = acceleration of gravity (default see gravity) 039 % rho = water density (default see wdensity) 040 % bet = 1, theta given in terms of directions toward which waves travel (default) 041 % -1, theta given in terms of directions from which waves come 042 % igam = 1, if z is measured positive upward from mean water level (default) 043 % 2, if z is measured positive downward from mean water level 044 % 3, if z is measured positive upward from sea floor 045 % thx = angle clockwise from true north to positive x-axis in degrees 046 % (default 90) 047 % thy = angle clockwise from true north to positive y-axis in degrees 048 % (default 0) 049 % kw = vector of wave numbers corresponding to angular frequencies, w 050 % (default calculated with w2k) 051 % 052 % Example: 053 % N=5000;dt=0.4;f0=0.1;th0=0;h=50; w0 = 2*pi*f0; 054 % t = linspace(0,1500,N)'; 055 % types = [1 4 5]; 056 % bfs = [1 0 0]; 057 % pos = zeros(3,3); 058 % eta0 = exp(-i*w0*t); 059 % [Hw Gwt] = tran(w0,th0,pos(1,:),'n',h); 060 % eta = real(Hw*Gwt*eta0)+0.001*rand(N,1); 061 % [Hw Gwt] = tran(w0,th0,pos(2,:),'n_x',h); 062 % n_x = real(Hw*Gwt*eta0)+0.001*rand(N,1); 063 % [Hw Gwt] = tran(w0,th0,pos(3,:),'n_y',h); 064 % n_y = real(Hw*Gwt*eta0)+0.001*rand(N,1); 065 % 066 % S = dat2dspec([t eta n_x n_y],[pos types',bfs'],h,256,101); 067 % wspecplot(S) 068 % 069 % See also dat2dspec, wdensity, gravity 070 071 % Reference 072 % Young I.R. (1994) 073 % "On the measurement of directional spectra", 074 % Applied Ocean Research, Vol 16, pp 283-294 075 076 % History 077 % revised pab jan2005 078 % -moved the code changing sensortypename to sensortypeid into a separate 079 % function sensortypeid.m 080 % revised pab 11.10.2002 081 % - fixed some bugs in transfer functions 082 % - reordered the transfer functions 083 % - 084 % revised pab 07.11.2001 085 % -Made a fix for w=0 where the transfer functions sometimes is singular. 086 % e.g., Hw = w.*ratio(zk,hk,-1,-1); return NaN for zk=0,hk=0 and w=0. 087 % This fix is strictly not correct for all cases, but suffices in most cases. 088 % -moved tran as a local function to a function in \wafo\multidim\ 089 % revised pab 14.06.2000 090 % - updated documentation 091 % - clearified directions 092 % - added def for surface curvature 093 % revised pab 06.01.2000 094 % - cleaned up documentation 095 % - added default values 096 % - speeded up computations 097 % based on transp by L. Borgman 098 099 % Default values 100 %~~~~~~~~~~~~~~~ 101 if nargin<2|isempty(theta), theta = 0; end 102 if nargin<3|isempty(pos), pos = [0 0 0]; end 103 if nargin<4|isempty(def), def = 1; end 104 if nargin<5|isempty(h), h = inf; end 105 if nargin<6|isempty(g), g = gravity;, end % accelleration of gravity 106 if nargin<7|isempty(rho), rho = wdensity; end % water density given in kg/m^3 107 if nargin<8|isempty(bet), bet = 1; end 108 if nargin<9|isempty(igam), igam = 1;end 109 if nargin<10|isempty(thx), thx = 90; end 110 if nargin<11|isempty(thy), thy = 0; end 111 if nargin<12|isempty(kw), kw = w2k(w,0,h); end % find the wave number as function of angular frequency 112 113 if ischar(def) 114 def1 = sensortypeid(def); 115 if isempty(def1)|length(def1)>1 116 error(sprintf('Invalid input def = %s',def)) 117 end 118 def = def1; 119 end 120 121 % make sure they have the correct orientation 122 theta = theta(:); 123 kw = kw(:).'; 124 w = w(:).'; 125 126 127 % Compute various fixed arrays 128 % ---------------------------- 129 Nt = length(theta); 130 Nf = length(w); 131 Oneth = ones(Nt,1); 132 Onef = ones(1,Nf); 133 134 135 136 % convert to from angle in degrees to radians 137 thxr = thx*pi/180; 138 thyr = thy*pi/180; 139 140 cthx = bet*cos(theta-thxr+pi/2); 141 %cthy = cos(theta-thyr-pi/2); 142 cthy = bet*sin(theta-thyr); 143 144 % Compute location complex exponential 145 % ------------------------------------ 146 147 x=pos(1);y=pos(2);z=pos(3); 148 %arg = bet*(x*cthx+y*cthy)*kw; 149 %arg = arg(:,Onef).*kw(Oneth,:); 150 151 % old call 152 %ee=cos(arg)-(sqrt(-1))*sin(arg); % exp(-i*k(w)*(x*cos(theta)+y*sin(theta)) size Nt X Nf 153 %new call 154 ee= exp((sqrt(-1)*(x*cthx+y*cthy))*kw); % exp(i*k(w)*(x*cos(theta)+y*sin(theta)) size Nt X Nf 155 156 157 if def > 8 158 % Set up z and h arguments 159 % ------------------------ 160 hk = kw*h; 161 switch igam 162 case 1, zk = kw*(h+z); % z measured positive upward from mean water level (default) 163 case 2, zk = kw*(h-z); % z measured positive downward from mean water level 164 otherwise,zk = kw*z; % z measured positive upward from sea floor 165 end 166 end 167 168 % Start computation of complex transformations 169 % --------------------------------------------- 170 switch def 171 case 1, % n = Eta = wave profile 172 Hw = Onef; 173 Gwt = ee; 174 175 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 176 % 177 % Vertical surface velocity and acceleration 178 % 179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 180 181 case 2, % n_t = Eta_t 182 Hw = w; 183 Gwt = -sqrt(-1)*ee; 184 case 3, % n_tt = Eta_tt 185 Hw = w.^2; 186 Gwt = -ee; 187 188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 189 % 190 % Surface slopes 191 % 192 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 193 case 4, % n_x = Eta_x = x-slope 194 Hw = kw; 195 Gwt = sqrt(-1)*cthx(:,Onef).*ee; 196 case 5, % n_y = Eta_y = y-slope 197 Hw = kw; 198 Gwt = sqrt(-1)*cthy(:,Onef).*ee; 199 200 201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 202 % 203 % Surface curvatures 204 % 205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 206 case 6, % n_xx = Eta_xx = Surface curvature (x-dir) 207 Hw = kw.^2; 208 Gwt = -cthx(:,Onef).^2.*ee; 209 case 7, % n_yy = Eta_yy = Surface curvature (y-dir) 210 Hw = kw.^2; 211 Gwt = -cthy(:,Onef).^2.*ee; 212 case 8, % n_xy = Eta_xy = Surface curvature (xy-dir) 213 Hw = kw.^2; 214 Gwt = -cthx(:,Onef).*cthy(:,Onef).*ee; 215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 216 % 217 % Pressure 218 % 219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 220 221 222 case 9, % pressure fluctuations 223 Hw = rho*g*ratio(zk,hk,1,1); % cosh(zk)/cosh(hk) 224 Gwt = ee; 225 226 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 227 % 228 % water particle velocities 229 % 230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 231 case 10, % U = x-velocity 232 Hw = w.*ratio(zk,hk,1,-1); % w*cosh(zk)/sinh(hk) 233 Gwt = cthx(:,Onef).*ee; % cos(theta)*ee 234 case 11, % V = y-velocity 235 Hw = (w.*ratio(zk,hk,1,-1)); % w*cosh(zk)/sinh(hk) 236 Gwt = cthy(:,Onef).*ee; % sin(theta)*ee 237 case 12, % W = z-velocity 238 Hw = w.*ratio(zk,hk,-1,-1); % w*sinh(zk)/sinh(hk) 239 Gwt = -sqrt(-1)*ee; %-? 240 241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 242 % 243 % water particle acceleration 244 % 245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 246 case 13 , % U_t = x-acceleration 247 Hw = (w.^2).*ratio(zk,hk,1,-1);% w^2*cosh(zk)/sinh(hk) 248 Gwt = -sqrt(-1)*cthx(:,Onef).*ee; %-? 249 case 14 , % V_t = y-acceleration 250 Hw = (w.^2).*ratio(zk,hk,1,-1);% w^2*cosh(zk)/sinh(hk) 251 Gwt = -sqrt(-1)*cthy(:,Onef).*ee; %-? 252 case 15, % W_t = z-acceleration 253 Hw = (w.^2).*ratio(zk,hk,-1,-1);% w*sinh(zk)/sinh(hk) 254 Gwt = -ee; 255 256 257 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 258 % 259 % water particle displacement 260 % 261 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 262 case 16, % X_p = x-displacement 263 Hw = ratio(zk,hk,1,-1); % cosh(zk)./sinh(hk) 264 Gwt = sqrt(-1)*cthx(:,Onef).*ee; 265 case 17, % Y_p = y-displacement 266 Hw = ratio(zk,hk,1,-1); % cosh(zk)./sinh(hk) 267 Gwt = sqrt(-1)*cthy(:,Onef).*ee; 268 case 18, % Z_p = z-displacement 269 Hw = ratio(zk,hk,-1,-1); % sinh(zk)./sinh(hk) 270 Gwt = ee; 271 otherwise 272 error('Unknown option 1<=def<=18') 273 end 274 275 % New call to avoid singularities. pab 07.11.200 276 % Set Hw to 0 for expressions w*ratio(z*k,h*k,1,-1)= 0*inf 277 Hw(isnan(Hw) | isinf(Hw)) = 0; 278 279 sgn = sign(Hw); 280 k0 = find(sgn<0); 281 if any(k0) % make sure Hw>=0 ie. transfer negative signs to Gwt 282 Gwt(:,k0) = -Gwt(:,k0); 283 Hw(k0) = -Hw(k0); 284 285 end 286 if igam==2, 287 %pab 09 Oct.2002: bug fix 288 % Changing igam by 2 should affect the directional result in the same way that changing eta by -eta! 289 Gwt = -Gwt; 290 end 291 292 293 if nargout>3 294 Hwt=Hw(Oneth,:).*Gwt; 295 end 296 %if Hw is wanted in size Nt x Nf uncomment this!!! 297 % Hw=Hw(Oneth,:); 298 299 return 300 301
Comments or corrections to the WAFO group