Home > wafo > multidim > tran.m

tran

PURPOSE ^

Computes transfer functions based on linear wave theory

SYNOPSIS ^

[Hw, Gwt, kw,Hwt,ee]=tran(w,theta,pos,def,h,g,rho,bet,igam,thx,thy,kw);

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003