function TF=orca_nba_xyz_run(BW,Tw,iimeld,iifft,iifig)
% 
% Call TF=orca_nba_ri_run(BW,Tw,iimeld,iifft,iifig);
%
% BW		= bandwidth about each center frequency (0==> default of 2*min(diff(f_cw)));
% Tw		= pos==>time window (1-2 s, e.g.), neg==>frequency bin spacing (default=2);
% iimeld	= meld NBA approximations together to form one transfer function across 
%			  multiple center frequencies (default=1);
% iifft		= 1 to output FFT file (default=0);
% iifig		= 1 to view output data structure TF (default=1);
%
%
%
% 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.
if nargin<1 | isempty(BW), BW=0; end
if nargin<2 | isempty(Tw), Tw=2; end
if nargin<3 | isempty(iimeld), iimeld=1; end
if nargin<4 | isempty(iifft), iifft=0; end
if nargin<5 | isempty(iifig), iifig=1; end
%
orca_global
%
cback_fig=Obox.figno;
%
type_opts=[1,0,0,0,0];
gen_opts=ORC.gen_opts;
iitl=0; iimf=1; iisig=0; iimt=0; iidc=1; iikn=1; iieig=0; iimf_dphi=0; iimf_im=1;
cw_opts=[iitl,iimf,iisig,iimt,iidc,iikn,iieig,0,0,0,0,0,0,iimf_dphi,iimf_im,0];
%
cf=adi_lim_eval(ORC.cw_freqs); ncf=length(cf);
if BW==0, 
	if ncf>1, BW=2*min(diff(cf)); else, disp('BW=0, but ncf=1.  Setting BW=30.'); BW=30; end
end
if Tw<0, Tw=-1/Tw; end
df=1/Tw; ndf=ceil(BW/df); ndf=floor(ndf/2)*2+1; ndfh=(ndf-1)/2; df_ax=[-ndfh:ndfh]*df;
%
zs=adi_lim_eval(ORC.zs_lims); nzs=length(zs);
zr=adi_lim_eval(ORC.zr_lims); nzr=length(zr);
rng_km=adi_lim_eval(ORC.rng_lims); rng=1000*rng_km; nrng=length(rng);
zmf_lims=unique([zs(:); zr(:)]);
iizs=zeros(1,nzs); iizr=zeros(1,nzr);
for jzs=1:nzs;
	[dzmin,iizs(jzs)]=min(abs(zs(jzs)-zmf_lims));
end
for jzr=1:nzr;
	[dzmin,iizr(jzr)]=min(abs(zr(jzr)-zmf_lims));
end
disp(iizs), disp(iizr),
%
iicomm=0;
% Output OPT file:
if strcmp(ORC.out_root,'auto'), out_root=[ORC.svp_root,'_NB']; else, out_root=ORC.out_root; end
cw_freqs=sort(cf); fmin=cw_freqs(1); fmax=cw_freqs(end);
orca_opt_write(out_root,iicomm,ORC.type_opts,ORC.gen_opts,cw_freqs,cw_opts(1:16), ...
		ORC.zs_lims,ORC.zr_lims,ORC.rng_lims,zmf_lims,ORC.bb_parms,ORC.bb_opts, ...
		ORC.var_parm1,ORC.var_parm2,ORC.kpl_parms,ORC.rc_parms);
% Run orca90:
svp_file=[ORC.adi_dir,ORC.svp_root,'.svp'];
opt_file=[out_root,'.opt'];
if 1,
	stat=orca90(svp_file,opt_file,out_root);
	if stat==0, disp('ORCA did not run succesfully'); return; else, disp('ORCA finished.'); end
end
%
[SD,status]=hdfsd_open([out_root,'.hdf']); if status<0, error('Reading SD'); end
jset_kn=find(strcmp(SD.data_name,'Re(kn)'));
if isempty(jset_kn), 
	error('Unexpected in orca_nba_ri_run: HDF file does not contain Re(kn), Im(kn).'); 
end
[KNR,status]=hdfsd_read(SD,[],jset_kn); if status<0, error('Reading KNR'); end
[KNI,status]=hdfsd_read(SD,[],jset_kn+1); if status<0, error('Reading KNI'); end
jset_vg=find(strncmp(SD.data_name,'Group Speed',11));
if isempty(jset_vg), 
	error('Unexpected in orca_nba_ri_run: HDF file does not contain Group Speed.'); 
end
[VG,status]=hdfsd_read(SD,[],jset_vg); if status<0, error('Reading VG'); end
%
jset_phi_re=find(strcmp(SD.data_name,'Re(phi)'));
jset_phi_im=find(strcmp(SD.data_name,'Im(phi)'));
if isempty(jset_phi_re) | isempty(jset_phi_im), error('Re/Im(phi) not found'); end
%
[PHI_re,status]=hdfsd_read(SD,[],jset_phi_re); if status<0, error('Reading PHI_re'); end
[PHI_im,status]=hdfsd_read(SD,[],jset_phi_im); if status<0, error('Reading PHI_im'); end
%
out_root=[out_root,'_NBA'];
TF.title=out_root;
TF.dlab='Transfer Function';
%
f0=cw_freqs; f0=round(f0/df)*df; 
nparm=size(KNR.data,3);
if iimeld,
	f1=f0(1) - ndfh*df; f2=f0(end) + ndfh*df;
	ndf_tot=round((f2-f1)/df) + 1;
	fmeld=linspace(f1,f2,ndf_tot);
	if nparm==1,
	    TF.data=zeros(nzr,nrng,ndf_tot,nzs);
	    TF.ndim=4;
	    TF.dims=[nzr,nrng,ndf_tot,nzs];
	    TF.labels={'Depth - m','Range - km','Frequency - Hz','Source Depth - m'};
	    TF.scales={zr,rng_km,fmeld,zs};
	else
	    TF.data=zeros(nzr,nrng,ndf_tot,nzs,nparm);
	    TF.ndim=5;
	    TF.dims=[nzr,nrng,ndf_tot,nzs,nparm];
	    TF.labels={'Depth - m','Range - km','Frequency - Hz','Source Depth - m','Parameter'};
	    TF.scales={zr,rng_km,fmeld,zs,parms};
	end
else
	if nparm==1,
	    TF.data=zeros(nzr,nrng,ndf,nzs,ncf);
	    TF.ndim=5;
	    TF.dims=[nzr,nrng,ndf,nzs,ncf];
	    TF.labels={'Depth - m','Range - km','Frequency - Hz','Source Depth - m','CF - Hz'};
	    TF.scales={zr,rng_km,df_ax,zs,f0};
	else
	    TF.data=zeros(nzr,nrng,ndf,nzs,ncf,nparm);
	    TF.ndim=6;
	    TF.dims=[nzr,nrng,ndf,nzs,ncf,nparm];
	    TF.labels={'Depth - m','Range - km','Frequency - Hz','Source Depth - m','CF - Hz','Parameter'};
	    TF.scales={zr,rng_km,df_ax,zs,f0,parms};
	end
end
TF.iidB=-1;
TF.vg_ref=zeros(ncf,nparm);
%

for jparm=1:nparm,
	for jcf=1:ncf,
        iim=find(~isnan(VG.data(jcf,:,jparm))); nmode=length(iim);
        [im_min,jm_ref]=min(KNI.data(jcf,iim,jparm)); vg_ref=VG.data(jcf,jm_ref,jparm);
        tau_norm=rng/vg_ref; TF.vg_ref(:,jparm)=vg_ref;
%
		phis=PHI_re.data(iizs,1:nmode,jcf,jparm);
		phir=PHI_re.data(iizr,1:nmode,jcf,jparm);
		vg=VG.data(jcf,1:nmode,jparm);
		kn=KNR.data(jcf,1:nmode,jparm) + i*KNI.data(jcf,1:nmode,jparm);
		f_ax=f0(jcf) + df_ax;
%
		disp(['Starting cf = ',int2str(jcf),' of ',int2str(ncf)]);
		if iimeld==0,
			TF.data(:,:,:,:,jcf,jparm)=orca_nba_ri(nmode,nrng,nzs,nzr,ndf,rng,f0(jcf),f_ax, ...
    			phis,phir,vg,kn,vg_ref);
		else
			f1=f0(jcf) - ndfh*df;
			jf1m=round((f1-fmeld(1))/df) + 1;
			f2=f0(jcf) + ndfh*df;
			jf2m=round((f2-fmeld(1))/df) + 1;
			if jcf==1,
				TF.data(:,:,jf1m:jf2m,:,jparm)=orca_nba_ri(nmode,nrng,nzs,nzr,ndf,rng,f0(jcf),f_ax, ...
    				phis,phir,vg,kn,vg_ref);
			else
				tf=orca_nba_ri(nmode,nrng,nzs,nzr,ndf,rng,f0(jcf),f_ax,phis,phir,vg,kn,vg_ref);
				jf3m=max(jf1m,jf2mx) + 1;
				jf1b=1; jf2b=round((f2x-f1)/df) + 1;
				nf_hann=jf2mx-jf1m+1;
				hann0=.5*(1 - cos((2*pi/(2*nf_hann))*[.5:1:nf_hann]));
				hann2=reshape(repmat(hann0,[nrng,1]),[1,nrng,nf_hann]);
				hann1=1 - hann2;
				for jzs=1:nzs,
					for jzr=1:nzr,
						TF.data(jzr,:,[jf1m:jf2mx],jzs,jparm)= ...
							hann1.*TF.data(jzr,:,[jf1m:jf2mx],jzs,jparm) + ...
							hann2.*tf(jzr,:,[jf1b:jf2b],jzs);
						TF.data(jzr,:,[jf3m:jf2m],jzs,jparm)=tf(jzr,:,[jf2b+1:ndf],jzs);
					end
				end
			end
			f2x=f2; jf2mx=jf2m;
		end
		
	end
end
%
str=['save ',out_root,' TF']; disp(str); eval(str);
%
if iifft,
	nrec=TF.dims(1); nfreq=TF.dims(3); nzp_f=2; ii_app=0;
	df=diff(TF.scales{3}(1:2)); f1=TF.scales{3}(1); f2=TF.scales{3}(end);
	FFT=fft_write1_open2(out_root,nrec,f1,f2,df,nzp_f,ii_app);
	jrec=[1:nrec];
	jday=0; t0=0; 
	jrec_parm=jrec;
	rng_m=TF.scales{2}*1000; c_ref=1500; w_ax=2*pi*TF.scales{3}.';
	for jseq=1:TF.dims(2),
		t_ref=rng_m(jseq)/c_ref - .05/df; exp_fac=exp(i*w_ax*t_ref);
		cspec=conj(reshape(TF.data(:,jseq,:,1,1),[nrec,nfreq]).');
		cspec=cspec.*repmat(exp_fac,[1,nrec]);
		count=fft_write1(FFT,cspec,jseq,jrec,jday,t0,jseq,jrec_parm);
	end
	fclose(FFT.fft_fid);
end

