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)

## CROSS-REFERENCE INFORMATION

This function calls:
 gravity returns the constant acceleration of gravity ratio compute ratio of hyperbolic functions sensortypeid Return ID for sensortype name w2k Translates from frequency to wave number wdensity Returns the water density error Display message and abort function. ischar True for character array (string). sprintf Write formatted data to string.
This function is called by:
 dat2dspec Estimates the directional wave spectrum from timeseries testmeasurements Creates a test case for measurement time series

## 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 %
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