function tf=orca_nba_xyz(xyz_src,xyz_rec,f0,f_ax,kn,vg,zmf,phi)
% 
% Call tf=orca_nba_xyz(xyz_src,xyz_rec,f0,f_ax,kn,vg,zmf,phi);
%
% tf		= [ndf,nrec,nsrc] complex transfer function
%
% xyz_src	= nsrc by 3 array of source positions in m;
% xyz_rec	= nrec by 3 array of receiver positions in m;
%
% Computes complex received spectrum tf(1:ndf,1:nrec,1:nsrc) using narrowband approx.
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
nrec=size(xyz_rec,1);
nsrc=size(xyz_src,1);
ndf=length(f_ax);
iim=find(~isnan(kn)); nmode=length(iim);
%
w_ax=2*pi*f_ax(:); w0=2*pi*f0;
fac=sqrt(2*pi)*exp(i*pi/4);
iim=[1:nmode];
tf=zeros(ndf,nrec,nsrc);
%
jzr=zeros(1,nrec);
for jrec=1:nrec,
	[dum,jzr(jrec)]=min(abs(zmf - xyz_rec(jrec,3)));
end
for jsrc=1:nsrc,
	xy_s=xyz_src(jsrc,1:2); 
	[dum,jzs]=min(abs(zmf - xyz_src(jsrc,3)));
	phi_src=phi(jzs,iim);
	for jrec=1:nrec,
		phi_rec=phi(jzr(jrec),iim);
		xy_r=xyz_rec(jrec,1:2);
		rng=sqrt(sum((xy_s - xy_r).^2));
		taux=rng./vg(iim);
		int_knr=kn(iim)*rng;
		Efac= exp( i * (int_knr - w0*taux) ) ./ sqrt(int_knr);	
% Do as matrices (ndf by nmode):
		Eshift=exp( (i*w_ax) * taux );
        phi_prod=fac*phi_src.*phi_rec;			% 1 by nmode
%%			size(phi_prod),size(Efac),size(Eshift),	keyboard
		tf_modes=repmat(phi_prod .* Efac,[ndf,1]) .* Eshift;
% Note: take conjugate before outputting to FFT file.  This convention agrees with
% phase of TL output by ORCA:
		tf(:,jrec,jsrc)=sum(tf_modes,2);
	end
end
