function [Hw, Gwt, kw,Hwt]=tran(w,theta,pos,def,h,g,rho,bet,igam,thx,thy,kw);
%TRAN computes matrix of transfer functions based on linear wave theory
%     of the system with input surface elevation, eta(x0,y0,t) = cos(kx*x0+ky*y0-w*t), 
%     and output Y determined by def and pos. 
%
% CALL:  [Hw Gwt kw Hwt] = tran(w,theta,pos,def,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])
%   def = Number defining the output 
%         1, wave profile                     (Eta)     (default)
%         2, x-slope of sea surface           (Eta_x)
%         3, y-slope of sea surface           (Eta_y)
%         4, x-velocity of water particle     (U)
%         5, y-velocity of water particle     (V)
%         6, z-velocity of water particle     (W)
%         7, x-acceleration of water particle (U_t)
%         8, y-acceleration of water particle (V_t)
%         9, z-acceleration of water particle (W_t)
%         10, pressure fluctuation about static mwl pressure  (P)
%         11, Surface curvature x-direction   (Eta_xx)
%         12, Surface curvature y-direction   (Eta_yy)
%         13, Surface curvature xy-direction  (Eta_xy)
%         14, Water particle displacement in x-direction
%         15, Water particle displacement in y-direction
%         16, Water particle displacement in z-direction
%         17, 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)
%
%
% See also: dat2dspec, wdensity, gravity

% Reference
% Young I.R. (1994)
% "On the measurement of directional spectra",
% Applied Ocean Research, Vol 16, pp 283-294

% History
% revised pab 07.11.2001
% -Made a fix for w=0 where the transfer functions sometimes is singular.
%  e.g., Hw  = w.*ratio(zk,hk,-1,-1); return NaN for zk=0,hk=0 and w=0.
%  This fix is strictly not correct for all cases, but suffices in most cases.
% -moved tran as a local function to a function in \wafo\multidim\
% revised pab 14.06.2000
% - updated documentation
% - clearified directions
% - added def for surface curvature
% revised pab 06.01.2000
%  - cleaned up documentation
%  - added default values
%  - speeded up computations
% based on transp by L. Borgman

% Default values
%~~~~~~~~~~~~~~~
if nargin<2|isempty(theta), theta = 0; end
if nargin<3|isempty(pos),   pos   = [0 0 0]; end
if nargin<4|isempty(def),   def   = 1; end  
if nargin<5|isempty(h),     h     = inf; end
if nargin<6|isempty(g),     g     = gravity;, end % accelleration of gravity
if nargin<7|isempty(rho),   rho   = wdensity; end % water density given in kg/m^3
if nargin<8|isempty(bet),   bet   = 1; end
if nargin<9|isempty(igam),  igam  = 1;end
if nargin<10|isempty(thx),  thx   = 90; end
if nargin<11|isempty(thy),  thy   = 0; end
if nargin<12|isempty(kw),   kw    = w2k(w,0,h); end % find the wave number as function of angular frequency


% make sure they have the correct orientation
theta = theta(:);
kw    = kw(:).';
w     = w(:).';


% Compute various fixed arrays
% ----------------------------
Nt    = length(theta);
Nf    = length(w);
Oneth = ones(Nt,1);
Onef  = ones(1,Nf);


% convert to from angle in degrees to radians
thxr = thx*pi/180;
thyr = thy*pi/180;

cthx = cos(theta-thxr+pi/2);
%cthy=cos(theta-thyr-pi/2);
cthy = sin(theta-thyr);

% Compute location complex exponential
% ------------------------------------

x=pos(1);y=pos(2);z=pos(3);
arg = bet*(x*cthx+y*cthy)*kw;
%arg = arg(:,Onef).*kw(Oneth,:);

% old call
%ee=cos(arg)-(sqrt(-1))*sin(arg); % exp(-i*k(w)*(x*cos(theta)+y*sin(theta)) size Nt X Nf
%new call
ee= exp((sqrt(-1))*arg);          % exp(i*k(w)*(x*cos(theta)+y*sin(theta)) size Nt X Nf


if def > 3
   % Set up z and h arguments
   % ------------------------ 
   hk=kw*h;
   switch igam
     case 1,   zk = kw*(h+z); % z measured positive upward from mean water level (default)
     case 2,   zk = kw*(h-z); % z measured positive downward from mean water level
     otherwise,zk = kw*z;     % z measured positive upward from sea floor
   end
end

% Start computation of complex transformations
% ---------------------------------------------
switch def
  case 1,  % wave profile
    Hw  = Onef;
    Gwt = ee;
  case 2,  % x-slope
    Hw  = kw;
    Gwt = sqrt(-1)*bet*cthx(:,Onef).*ee;
  case 3,  % y-slope
    Hw  = kw;
    Gwt = sqrt(-1)*bet*cthy(:,Onef).*ee;
  case 4,  % x-velocity
    Hw  = w.*ratio(zk,hk,1,-1);
    Gwt = bet*cthx(:,Onef).*ee; % cos(theta)*ee 
  case 5,  % y-velocity
    Hw  = (w.*ratio(zk,hk,1,-1));
    Gwt = bet*cthy(:,Onef).*ee; 
  case 6,  % z-velocity
    Hw  = w.*ratio(zk,hk,-1,-1); 
    Gwt = -sqrt(-1)*ee;                 
  case 7 , % x-acceleration
    Hw  = (w.^2).*ratio(zk,hk,1,-1);
    Gwt = -sqrt(-1)*bet*cthx(:,Onef).*ee;
  case 8  ,% y-acceleration
    Hw  = (w.^2).*ratio(zk,hk,1,-1);
    Gwt = -sqrt(-1)*bet*cthy(:,Onef).*ee;
  case 9,  % z-acceleration
    Hw  =(w.^2).*ratio(zk,hk,-1,-1);
    Gwt =-ee;
  case 10, % pressure fluctuations
    Hw  = rho*g*ratio(zk,hk,1,1);
    Gwt = ee;
  case 11, % Surface curvature (x-dir)
    Hw = kw.^2;
    Gwt = -bet*cthx(:,Onef).^2.*ee; % bet???
  case 12, % Surface curvature (y-dir)
    Hw = kw.^2;
    Gwt = -bet*cthy(:,Onef).^2.*ee; % bet???
  case 13, % Surface curvature (xy-dir)
    Hw = kw.^2;
    Gwt = -bet*cthx(:,Onef).*cthy(:,Onef).*ee; % bet???
  case 14, % x-displacement
    Hw  = ratio(zk,hk,1,-1);  %cosh(zk)./sinh(hk) 
    Gwt = -sqrt(-1)*bet*cthy(:,Onef).*ee;
  case 15, % y-displacement
    Hw  = ratio(zk,hk,-1,-1); % sinh(zk)./sinh(hk)
    Gwt = -sqrt(-1)*bet*cthx(:,Onef).*ee;
  case 16, % z-displacement
    Hw  = rho*g*ratio(zk,hk,1,1); % rho*g*cosh(zk)./cosh(hk) 
    Gwt = ee;
  case 17, % dummy
    Hw  = zeros(1,Nf);
    Gwt = zeros(Nt,Nf);
  otherwise
    error('Unknown option 1<=def<=17')
 end
 
 % New call to avoid singularities. pab 07.11.200
 % Set Hw to 0 for expressions w*ratio(z*k,h*k,1,-1)= 0*inf
Hw(isnan(Hw) | isinf(Hw)) = 0;

sgn = sign(Hw);
if any(sgn<0) % make sure Hw>=0 ie. transfer negative signs to Gwt
  Gwt = sgn(Oneth,:).*Gwt;
  Hw  = abs(Hw);
end


if nargout>3
 Hwt=Hw(Oneth,:).*Gwt;
end
%if Hw is wanted in size Nt x Nf uncomment this!!!
%   Hw=Hw(Oneth,:); 

return