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) 

%Edits were made on 3/10/09 to fix a bug in the cutoff limits in the KZ
%and KZi parts of the transfer function.  NOTE: there is a cutoff set of .4
%for both the real and imaginary parts, but this ONLY applies to high
%frequencies.  A cutoff is in place for the imaginary KZi at .05Hz to limit
%the effects of low frequency noise found in some data sets (at .01-.05Hz).
%This cutoff is set at a designated frequency to avoid losing variance at
%frequencies above .05Hz (20 second waves).



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)


%what is set for the cutoff value to be used at high frequencies
cutoff=.4;

%THE REAL PART
KZ=cosh(wns*z)./sinh(wns*depth);
FLKZ=ffreqs.*KZ;

%find the cutoff and replace the values below it
cutindex=find(FLKZ < cutoff);
FLKZ(cutindex)=cutoff;


%THE IMAGINARY PART
KZi=sinh(wns*z)./sinh(wns*depth);
FLKZi=ffreqs.*KZi;

%compute the max index
[maxvalue,maxindex]=max(FLKZi);
%find the cuttoff and replace the values for high frequencies
cutindex=find(FLKZi(maxindex:end) < cutoff);
cutindex=cutindex+(maxindex-1);
FLKZi(cutindex)=cutoff;
%now insert cutoff values at frequencies < .06Hz
lowindex=find(ffreqs < (.06*2*pi));
lowvalue=FLKZi(max(lowindex)+1);
FLKZi(1:max(lowindex))=lowvalue;


%now put them together to create the transfer matrix
trm=(FLKZ*(sin(a)*cos(dirs-B))-(1i*FLKZi)*A1);