function tf=orca_nba_ri(nmode,nrng,nzs,nzr,ndf,rng,f0,f_ax, ...
    phis,phir,vg,kn,vg_ref)
%
% Computes complex received spectrum tf(1:ndf,1:nzr,1:nzs) 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.
w_ax=2*pi*f_ax(:); w0=2*pi*f0;
fac=sqrt(2*pi)*exp(i*pi/4);
iim=[1:nmode];
tf=zeros(nzr,nrng,ndf,nzs);
%
for jrng=1:nrng,
% 1 by nmode:
%%	taux=rng(jrng)./vg(iim) - rng(jrng)/vg_ref;
% Don't subtract out reference travel time here:
	taux=rng(jrng)./vg(iim);
	int_knr=kn(iim)*rng(jrng);
	Efac= exp( i * (int_knr - w0*taux) ) ./ sqrt(int_knr);	
% Do as matrices (ndf by nmode):
	Eshift=exp( (i*w_ax) * taux );
	for jzs=1:nzs,
		for jzr=1:nzr,
	        phi_prod=fac*phis(jzs,iim).*phir(jzr,iim);			% 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(jzr,jrng,:,jzs)=sum(tf_modes,2);
		end
	end
end

