root/adcp/trunk/adcp/diwasp_1_1GD/private/radial.m

function trm=radial(ffreqs,dirs,wns,z,depth,xpos,ypos);

%This transfer function takes the along beam radial velocities from an ADCP
%and transfers them to surface displacement for use in DIWASP.
%This was added to DIWASP on 6/22/07
% Modified from Hoshimoto (1997) 



transH=.4; % transducer face height off the bottom
dist=(xpos^2+ypos^2)^.5; %horizontal distance from origin
a=atan(dist/(z-transH));%beam angle from the vertical

% compute the axis angle
if xpos > 0 & ypos > 0 % 1st quadrant
    B= acos(xpos/dist);
elseif xpos < 0 & ypos < 0 % 3 quadrant
    B= pi + acos(abs(xpos/dist));
elseif xpos < 0 & ypos > 0 % 2nd quadrant 
    B= pi/2 + asin(abs(xpos/dist));
elseif xpos > 0 & ypos < 0 % 4th quadrant
    B = 3*pi/2 + asin(abs(xpos/dist));
end

%note: Dirs and B are both measured from the x-axis with positive angles
%being counterclockwise

A1=(ones(1,length(dirs))*cos(a));%a vector of length dirs of the cos(a)

FL=ffreqs*ones(1,length(dirs));%Create a vector of length ffreqs/length dirs so that we can multiply
%through in the trm function 

KZ=cosh(wns*z)./sinh(wns*depth);
%include a maximum cuttoff for the velocity response function
KZ(find(KZ<0.1))=0.1;
KZ(find(isnan(KZ)))=1;

KZi=sinh(wns*z)./sinh(wns*depth);
%include a maximum cuttoff for the velocity response function
KZi(find(KZi<0.1))=0.1;
KZi(find(isnan(KZi)))=1;


trm=FL.*(KZ*(sin(a)*cos(dirs-B))-(i*KZi)*A1);