function [RGRAM,FFT]=rgram_3d(fft_root,REP,T_lims,jrec,f_lims,tau_lims,parms, ...
	jrec_rep,jseq_rep,iifig,clim);
%
% Call [RGRAM,FFT]=rgram_3d(fft_root,REP,T_lims,jrec,f_lims,tau_lims,parms, ...
%	jrec_rep,jseq_rep,iifig,clim);
%
% fft_root	= Root of FFT file (suffix = _fft) ([] to be prompted);
% REP	= Root of FFT file to crosscorrelate against;
% T_lims		= FFT sequences to output ([] for all);  Column vector of sequences, or
%			  [Tmin,Tmax], where T in hod for iihr=1, or minutes into FFT file 
%			  for iihr=0, or range in km for iihr=3;
% jrec		= List of receiver numbers to output;
% f_lims	= [fmin,fmax], or column vector of frequency indices, to use in IFFT;
% tau_lims	= [taumin,taumax] time delay interval to output ([] for all) in s;
% parms		= [dfnorm,iihr,iienv,fzpad,iihann,iinorm,ns_avg], where:  
% 	dfnorm	= 1 for data FFT, 0 for simulation, where range is used on vertical axis;
% 	iihr	= vertical axis units: 0 = min into FFT, 1 = HOD, 2 = range in km;
%	iienv	= flag to envelope (0=no;1=|env|;2=dB(|env|);3=arg(env)/pi;4=unwrap(arg(env))/pi;
%			  -1=complex env);
%	fzpad 	> 0 to zero pad by 2^fzpad (fzpad<0 ==> reduce fs),
%			< 0 to reduce sample rate and FFT size (when enveloping);
% 	iihann	= 1 to put a Hann window in the frequency domain before inverse FFTing (def=1);
%	iinorm	= 1 to apply correlation norm to replica correlation;
%	ns_avg	= number of sequences to average before outputing (default = 1);
% jrec_rep	= List of receivers (shifts) in REP FFT to use ([]=all);
% jseq_rep	= List of sequences (chirps) in REP FFT to use ([]=all);
% iifig		= 1 to call viewn(RGRAM);
% clim		= [cmin,cmax] for colorbar.
%
% OUTPUT is structure RGRAM.  Use viewn to produce image or line plots.
% RGRAM.scales=		= {tau, sequence or time, receiver/beam, replica FFT receiver, replica FFT seq};
% RGRAM.data	= array(1:ntau,1:nseq,1:nrec,1:nrec_rep,1:nseq_rep);
%
%
% 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, fft_root=[]; end
if nargin<2, REP=[]; end
if nargin<3, T_lims=[]; end
if nargin<4, jrec=[]; end
if nargin<5, f_lims=[]; end
if nargin<6, tau_lims=[]; end
if nargin<7 | isempty(parms), parms=[1 1 2 -3 1 0 1]; end
if nargin<8, jrec_rep=[]; end
if nargin<9, jseq_rep=[]; end
if nargin<10, iifig=1; end
if nargin<11, clim=[]; end
dfnorm=parms(1); iihr=parms(2); iienv=parms(3); fzpad=parms(4); 
if length(parms)>=5, iihann=parms(5); else, iihann=1; end
if length(parms)>=6, iinorm=parms(6); else, iinorm=0; end
if length(parms)>=7, ns_avg=parms(7); else, ns_avg=1; end
%
fft_file=[fft_root,'_fft'];
if exist(fft_file)~=2,
	disp(['Note: ',fft_file,' not found locally.']);
	[fft_file,path_name]=uigetfile('*_fft','Please Choose _fft File');
	if fft_file==0, RGRAM=[]; return; end
	lfft=length(fft_file); fft_root=fft_file(1:lfft-4);
	if ~strcmp(path_name,pwd), eval(['cd ',path_name]); end
end
%
% Open data FFT file:
FFT=fft_read1_open(fft_root,dfnorm);
%
% Open replica correlation FFT file:
%%REP=fft_read1_open(rep_root,0);
% Compute number of replicas to use:
if isempty(jrec_rep), jrec_rep=[1:REP.nrec]; end
if isempty(jseq_rep), jseq_rep=[1:REP.nseq]; end
nrec_rep=length(jrec_rep); nseq_rep=length(jseq_rep); nrep=nrec_rep*nseq_rep;
%
if abs(FFT.df-REP.df)/FFT.df > .001,
	error(['USER ERROR: df of FFT file and replica file do not match: ', ...
		num2str(FFT.df),',',num2str(REP.df)]);
end
%
if isempty(f_lims), f_lims=[FFT.f1,FFT.f2]; end
if any(size(f_lims)~=[1 2]),
	disp('Bad frequency specification: should be [fmin,fmax] or []'); return;
end
% Find frequency overlap between data and replica FFT files:
f_ovl=[max([FFT.f1,REP.f1,f_lims(1)]), min([FFT.f2,REP.f2,f_lims(2)])];
if f_ovl(2) <= f_ovl(1), 
	disp('No overlap between data FFT, replica FFT, and desired frequency bands FFTs'); 
	disp([FFT.f1,FFT.f2,REP.f1,REP.f2,f_lims]);
	return; 
end
%
fax=(FFT.imap(:,1) - 1) * FFT.df;
jfreq=find(fax>=f_ovl(1)-FFT.df/2 & fax<=f_ovl(2)+FFT.df/2);
nfreq=length(jfreq); fax_get=fax(jfreq); 
%
fax_rep=(REP.imap(:,1) - 1) * REP.df;
jfreq_rep=find(fax_rep>=f_ovl(1)-REP.df/2 & fax_rep<=f_ovl(2)+REP.df/2);
nfreq_rep=length(jfreq_rep);
%
if isempty(jrec), jrec=[1:FFT.nrec]; end
%
fac=2^fzpad;
dt=fac*FFT.dt;
%
[nr,nc]=size(T_lims);
if isempty(T_lims), 
	jseq=[1:FFT.nseq];
elseif nr==1 & nc==2,
% convert T_lims in min after start of FFt file to hod;
	if iihr==0, T_lims=mod(FFT.moy0+T_lims,1440)/60; end
	if iihr==2,
		jseq=find(FFT.range>=T_lims(1) & FFT.range<=T_lims(2));
	else	
		jseq=find(FFT.hodax>=T_lims(1) & FFT.hodax<=T_lims(2));
	end
else,
	jseq=T_lims(:)';
end
if isempty(jseq),
	T_lims,FFT.hodax(1),FFT.hodax(length(FFT.hodax)),
	disp('USER ERROR: No FFT sequences to be processed.'); keyboard; return;
end
%	
nrec=length(jrec); nseq=length(jseq); jseq_out=jseq(ns_avg:ns_avg:nseq);
nseq_out=floor(nseq/ns_avg);
RGRAM.fft_root=fft_root;
RGRAM.rep_root=REP.fft_root;
RGRAM.title=[fft_root,'*',REP.fft_root];
RGRAM.ndim=5; 
RGRAM.dims=zeros(1,5);
RGRAM.scales=cell(1,5);
%
if ns_avg<1 | round(ns_avg)~=ns_avg,
	disp(['USER ERROR: ns_avg must be integer>=1: ',int2str(ns_avg)]);
	return
end
%
RGRAM.data=[];
RGRAM.labels={'Delay Time - s', 'Time - min', 'Receiver/Beam', 'Replica Rec', 'Replica Seq'};
if iihr==2,
	RGRAM.labels{2}='Range - km'
elseif iihr==1,
	RGRAM.labels{2}=['Time - hour of day ',int2str(FFT.jday0)];
end
RGRAM.dlab='Pressure - \muPa';
if iinorm, RGRAM.dlab='Correlation Coefficient'; end
%
if iihann,
	delf=fax_get(nfreq) - fax_get(1);
	hann=.5*(1 - cos((2*pi/delf)*(fax_get-fax_get(1)))); 
% Hanning window widens the mainlobe, reduces it by 1.76 dB, and lowers the sidelobes.
	hann_norm=sqrt(sum(hann.^2)/length(hann));
% hann_arr is nfreq by nrec:
	hann_arr=repmat((hann/hann_norm),[1,nrec]);
end
%
rspec=zeros(nfreq,nrec_rep,nseq_rep); rep_spec=zeros(nfreq,nrec_rep);
% Read in Replica correlation FFTs:
for js_rep=1:nseq_rep,
	js_repx=jseq_rep(js_rep);
	[XHr,rep_spec(:,:)]=fft_read1(REP,js_repx,jrec_rep,jfreq_rep);
	rep_seq_ax(js_rep)=XHr.t0;
% Apply correlation norm always, then norm to signal strength for iinorm=0:
	rnorm=sqrt(sum(zmagsq(rep_spec),1));
	rep_spec(:,:)=rep_spec./repmat(rnorm,[nfreq,1]);
	rspec(:,:,js_rep)=rep_spec;
end
%num2str([rspec(1:10,1,1)],8)
%num2str([rspec(nfreq-10:nfreq,1,1)],8)
%nfreq,
RGRAM.scales{4}=XHr.recs;
RGRAM.scales{5}=rep_seq_ax;
%
fac=2^fzpad;
nfftx=FFT.nfft*fac; fsx=FFT.fs*fac; nT_intx=fac*FFT.nT_int;
if isempty(tau_lims), tau_lims=[0,FFT.Tfft-1/fsx]; end
tau_lims(1)=max(tau_lims(1),0); tau_lims(2)=min(tau_lims(2),FFT.Tfft-1/fsx);
jt=round(tau_lims*fsx)+1; iit=jt(1):jt(2); ntau=length(iit);
t=[iit-1]'*(1/fsx); RGRAM.scales{1}=t;
% In case either of tau_lims < 0:
iineg=find(iit<1); iit(iineg)=nfftx + iit(iineg);
% Initialize inverse FFT array and mapping indices:
fft_full=zeros(nfftx,nrec) + i*zeros(nfftx,nrec);
if iienv==0,
	imap1=FFT.imap(jfreq,1); imap2=nfftx + 2 - imap1;
	pser=zeros(nfftx,nrec); x=zeros(ntau,nrec);
else
% Shift band so f_bb is centered on f=0:
	f_bb=(REP.f1 + REP.f2)/2;
	disp('f_bb from REP.f1,REP.f2: '); disp(f_bb); 
	jf_bb=round((f_bb - fax_get(1))/FFT.df) + 1,
	if jf_bb<=0 | jf_bb>length(fax_get),
		error(['Replica FFT file not in same band as data FFT file.  Aborting.']); return
	end
	if abs(f_bb-fax_get(jf_bb)) > .001,
		disp(['NOTE: basebanding bin not equal to desired f_bb: ', ...
		num2str([f_bb,fax_get(jf_bb),FFT.df])])
	end
	pser=fft_full; x=zeros(ntau,nrec) + i*ones(ntau,nrec);
	iipos=[jf_bb:nfreq]; npos=length(iipos); iipos_bb=[1:npos];
	iineg=[1:jf_bb-1]; nneg=length(iineg); iineg_bb=[nfftx-nneg+1:nfftx];
end
%
% Correlator normalization: Use same expression as in cgram_3d, but correct
% for possibility that one time series may be heavily zero-padded.
% (This has been checked out on replica FFTs with single and multiple PRIsl)
max_fac=max(FFT.nT_int/FFT.nfft,REP.nT_int/REP.nfft);
min_fac=min(FFT.nT_int/FFT.nfft,REP.nT_int/REP.nfft);
% nfftx term cancels (1/nfftx) in Matlab inverse FFT.  Factor of 2 accounts for fact
% that the two spectra are being normed by only their positive frequency halves,
% which results in two missing 1/sqrt(2) factors:
fft_norm=(nfftx/2)*sqrt(max_fac/min_fac);
% In case of enveloping, account for negative frequencies that are filtered out:
% This term IS needed here to make correlations unity when enveloping:
if iienv, fft_norm=2*fft_norm; end
%
% When the replica is a perfect match to at least part of the data in rgram_3d
% and the data are passed through sgram_3d to get SGRAM)
% AND iinorm is set to 0:
% zmax(RGRAM.data,RGRAM.iidB) = zrms(SGRAM.data,SGRAM.iidB)
%
% When the correlator norm is applied (iinorm=1) for the above case,
% zmax(RGRAM.data,RGRAM.iidB) = 0 dB (unity correlation)
%
ns_out=0; 
for js=1:nseq,
	if js==1,
        RGRAM.dims=[ntau,nseq_out,nrec,nrec_rep,nseq_rep];
		RGRAM.data=zeros(ntau,nseq_out,nrec,nrec_rep,nseq_rep);
		if ns_avg>1, cspec_avg=zeros(nfreq,nrec,nrec_rep,nseq_rep); end
		cspec=zeros(nfreq,nrec); cspec0=cspec;
	end
	js_out=jseq(js);
	[XH,cspec(:,:)]=fft_read1(FFT,js_out,jrec,jfreq);
	if js==1, RGRAM.scales{3}=XH.recs; end
% Must hann right away; Otherwise, correlations won't be normalized:
	if iihann, cspec=cspec.*hann_arr; end
% Apply correlator norm to cspec if desired:
	if iinorm,
		cnorm=sqrt(sum(zmagsq(cspec),1));
		cspec(:,:)=cspec./repmat(cnorm,[nfreq,1]);
	end
%
	cspec0(:,:)=cspec(:,:);
% Loop over replica FFTs:
	for jseq_rep=1:nseq_rep,
		for jrec_rep=1:nrec_rep,
			jrep=(jseq_rep-1)*nrec_rep + jrec_rep;
%%			jjrep=[jrep:nrep:nrec_out]; 		% nrec long;
% Take conj of rspec so that time delay is tau(sspec) - tau(rspec), so
% if cspec has longer travel time and rspec has min travel time, peak
% will be at positive tau on output:
%%	data_intensity=10*log10(sum(zmagsq(cspec),1)),
%%	rep_intensity=10*log10(sum(zmagsq(rspec),1)),
			for jr=1:nrec,
% cspec is nfreq by nrec:
				cspec(:,jr)=cspec(:,jr).*conj(rspec(:,jrec_rep,jseq_rep));
			end
%%	mf_intensity=10*log10(sum(zmagsq(cspec),1)),
% Average over sequence number (meant for averaging cross_spectra):
			if ns_avg>1,
				cspec_avg(:,:,jrec_rep,jseq_rep)= ...
				cspec_avg(:,:,jrec_rep,jseq_rep) + cspec;
				if mod(js,ns_avg)==0,
					cspec(:,:)=cspec_avg(:,:,jrec_rep,jseq_rep)/ns_avg;
					cspec_avg(:,:,jrec_rep,jseq_rep)=0;
				end
			end
			if mod(js,ns_avg)==0,
				if jrep==1, ns_out=ns_out + 1; end
% Do inverse FFTs:
				if iienv==0,
					fft_full(imap1,:)=fft_norm*cspec;
					fft_full(imap2,:)=conj(fft_full(imap1,:));
					pser(:,:)=ifft(fft_full,nfftx,1);	% nfftx by nrec;
					x(:,:)=real(pser(iit,:));			% ntau by nrec;
				else
					fft_full(iipos_bb,:)=fft_norm*cspec(iipos,:);
					fft_full(iineg_bb,:)=fft_norm*cspec(iineg,:);
					pser(:,:)=ifft(fft_full,nfftx,1);	% nfftx by nrec;
					x(:,:)=pser(iit,:);					% ntau by nrec;
				end
				if iienv==0,
%%					x(:,:)=real(x);		% Already done above;
				elseif iienv==1,
					x(:,:)=abs(x);
				elseif iienv==2,
					x(:,:)=10*log10(max(1e-30,zmagsq(x)));
				elseif iienv==3,
					x(:,:)=angle(x)/pi;
				elseif iienv==4,
					x(:,:)=unwrap(angle(x))/pi;
				end
				RGRAM.data(:,ns_out,:,jrec_rep,jseq_rep)=x(:,:);
				cspec(:,:)=cspec0;
			end
		end
%%		if ~isdata, RGRAM.t_ref(ns_out)=XH.t0; end
		if iihr==1,
			xax(js)=XH.hod;
		elseif iihr==0,
			xax(js)=XH.moy + XH.t0/60 - FFT.moy0;
		end
	end
%
%%	[t,pser]=myffti(cspec,FFT,tau_lims,0,fzpad,iienv);
	if mod(js,10)==1, disp(['Finished sequence ',int2str(js),' of ',int2str(nseq)]); end
end
%
RGRAM.scales{2}=xax;
%%nseq_out,nrec,nrep,
%%keyboard
%
if isempty(clim),
	cmax=max(RGRAM.data(:));
	if iienv==0
		cmin=min(RGRAM.data(:));
	elseif iienv==1,
		cmin=0;
	elseif iienv==2,
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(RGRAM.data(:))/5)*5,cmax-60);
		RGRAM.dlab='Pressure - dB//\muPa^2';
	elseif iienv==3,
		cmin=-1; cmax=1;
		RGRAM.dlab='arg(env)/pi';
	elseif iienv==4,
		cmin=min(RGRAM.data(:));
		RGRAM.dlab='unwrap[arg(env)]/pi';
	elseif iienv==-1,
		cmax=max(dB(RGRAM.data(:)));
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(dB(RGRAM.data(:)))/5)*5,cmax-60);
		RGRAM.dlab='Complex pressure - \muPa';
	end
	clim=[cmin,cmax];
elseif iinorm,
	if iienv==2,
		clim=[clim(1)-clim(2),0];
	end
end
%
RGRAM.dt=RGRAM.scales{1}(min(ntau,2))-RGRAM.scales{1}(1);
RGRAM.pri=REP.pri; %%RGRAM.nT_int=round(REP.nT_int*REP.dt/RGRAM.dt);
RGRAM.T_rep=REP.nT_int*REP.dt;
RGRAM.T_lims=T_lims; RGRAM.tau_lims=tau_lims;
RGRAM.Tfft=FFT.Tfft;
RGRAM.T_step=FFT.T_step; RGRAM.fft_hann=FFT.iihann;
RGRAM.fzpad=fzpad; RGRAM.iienv=iienv;
RGRAM.clim=clim;
if iienv==2, RGRAM.iidB=1; elseif iienv==-1, RGRAM.iidB=-1; else, RGRAM.iidB=0; end
%
fclose(FFT.fft_fid);
%
if iifig & ~isempty(RGRAM.data),
	viewn(RGRAM,[4 1 2]); 
end
%
%%%
%
