function [CGRAM,FFT]=cgram_3d(fft_root1,fft_root2,Tlims,jrec,jfreq,tlims,parms, ...
	iifig,clim);
%
% Call [CGRAM,FFT]=cgram_3d(fft_root,fft_root2,Tlims,jrec,jfreq,tlims,parms, ...
%	iifig,clim);
%
% fft_root1	= Root of FFT file (suffix = _fft) ([] to be prompted);
% fft_root2	= Root of second FFT file (suffix = _fft) ([] for same as fft_root1);
% Tlims		= 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 crosscorrelate: 
%			  [jrec1_fft1 jrec1_fft2; ... jrecn_fft1 jrecn_fft2];
% jfreq		= [fmin,fmax], or column vector of frequency indices, to use in IFFT;
% tlims		= [tmin,tmax] time interval to output ([] for all) in seconds;
% parms		= [T_exp,iiihr,iienv,fzpad,iihann,ns_avg], where:  
% 	T_exp	= Exponential time factor in s for spectral whitening (0=no whitening, -1=phase transform);
% 	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);
%	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);
% iifig		= 1 to call view5(CGRAM);
% clim		= [cmin,cmax] for colorbar.
%
% OUTPUT is structure CGRAM.  Use view5 to produce image or line plots.
% CGRAM.dims		= [time, sequence or time, correlation pair];
%
if nargin<1, fft_root1=[]; end
if nargin<2, fft_root2=[]; end
if nargin<3, Tlims=[]; end
if nargin<4, jrec=[]; end
if nargin<5, jfreq=[]; end
if nargin<6, tlims=[]; end
if nargin<7 | isempty(parms), parms=[1 1 2 -3 0]; end
if nargin<8, iifig=1; end
if nargin<9, clim=[]; end
T_exp=parms(1); iihr=parms(2); iienv=parms(3); fzpad=parms(4); 
if length(parms)>=5, iihann=parms(5); else, iihann=0; end
if length(parms)>=6, ns_avg=parms(6); else, ns_avg=1; end
%
fft_file1=[fft_root1,'_fft'];
if exist(fft_file1)~=2,
	disp(['Note: ',fft_file1,' not found locally.']);
	[fft_file1,path_name]=uigetfile('*_fft','Please Choose _fft File');
	if fft_file1==0, CGRAM=[]; return; end
	lfft1=length(fft_file1); fft_root1=fft_file1(1:lfft1-4);
	eval(['cd ',path_name]);
end
%
FFT=fft_read1_open(fft_root1,1);
if isempty(fft_root2),
	FFT2=FFT;
else
	FFT2=fft_read1_open(fft_root2,1);
end
fax=(FFT.imap(:,1) - 1) * FFT.df;
[nr,nc]=size(jfreq);
if nr==1 & nc==2,
	jfreq=find(fax>=jfreq(1)-FFT.df/2 & fax<=jfreq(2)+FFT.df/2);
elseif isempty(jfreq), 
	jfreq=[1:FFT.nfreq];
end
nfreq=length(jfreq); fax_get=fax(jfreq); 
if isempty(jrec), 
	jrec1=[1:FFT.nrec]'; jrec2=[1:FFT2.nrec]';
	if FFT.nrec==FFT2.nrec | FFT2.nrec==1,
		jrec=repmat([1:FFT.nrec]',[1,2]);
	elseif FFT.nrec==1,
		jrec=[[1:FFT2.nrec]',ones(FFT2.nrec,1)];
	elseif FFT2.nrec==1,
		jrec=[[1:FFT.nrec]',ones(FFT.nrec,1)];
	end
end
%
fac=2^fzpad;
dt=fac*FFT.dt;
%
[nr,nc]=size(Tlims);
if isempty(Tlims), 
	jseq=[1:FFT.nseq];
elseif nr==1 & nc==2,
% convert Tlims in min after start of FFt file to hod;
	if iihr==0, Tlims=mod(FFT.moy0+Tlims,1440)/60; end
	if iihr==2,
		jseq=find(FFT.range>=Tlims(1) & FFT.range<=Tlims(2));
	else	
		jseq=find(FFT.hodax>=Tlims(1) & FFT.hodax<=Tlims(2));
	end
else,
	jseq=Tlims(:)';
end
if isempty(jseq),
	Tlims,FFT.hodax(1),FFT.hodax(length(FFT.hodax)),
	disp('USER ERROR: No FFT sequences to be processed.'); keyboard; return;
end
%
if ns_avg<1 | round(ns_avg)~=ns_avg,
	disp(['USER ERROR: ns_avg must be integer>=1: ',int2str(ns_avg)]);
	return
end
%
nrec=size(jrec,1); nseq=length(jseq); jseq_out=jseq(ns_avg:ns_avg:nseq);
nseq_out=floor(nseq/ns_avg);
%
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
%
cnorm1=zeros(1,nrec); cnorm2=zeros(1,nrec); cnorm=zeros(nfreq,nrec);
%
fac=2^fzpad,
nfftx=FFT.nfft*fac, fsx=FFT.fs*fac; nT_intx=fac*FFT.nT_int
if isempty(tlims), tlims=[-FFT.Tfft/2,FFT.Tfft/2-1/fsx]; end
tlims(1)=max(tlims(1),-FFT.Tfft/2); tlims(2)=min(tlims(2),FFT.Tfft/2-1/fsx);
jt=round(tlims*fsx)+1; 
iit0=[jt(1):jt(2)]; ntime=length(iit0);
iit=iit0;
%%iit=mod(iit0-1,nfftx)+1; min(iit), max(iit),
t=[iit0-1]'*(1/fsx); CGRAM.y=t;
% In case either of tlims < 0:
iineg=find(iit<1); iit(iineg)=nfftx + iit(iineg);
fft_full=zeros(nfftx,nrec) + i*zeros(nfftx,nrec);
if iienv,
% Shift band so it is centered on f=0:
	iih=ceil(nfreq/2); iix=nfftx-nfreq+iih + 1;
	pser=fft_full; x=zeros(ntime,nrec) + i*ones(ntime,nrec);
else
	imap1=FFT.imap(jfreq,1); imap2=nfftx + 2 - imap1;
	if any(imap2<1), error('Bad imap2: probably illegal ZP factor'); end
	pser=zeros(nfftx,nrec); x=zeros(ntime,nrec);
end
%
CGRAM.fft_root1=fft_root1;
CGRAM.fft_root2=fft_root2;
CGRAM.title=[fft_root1,' x ',fft_root2];
CGRAM.jrec=jrec;
%
CGRAM.ndim=3;
CGRAM.dims=[ntime,nrec,nseq_out];
CGRAM.scales={t, [1:nrec], FFT.hodax(jseq_out)};
CGRAM.labels={'Delay Time - s', 'Receiver Pair', ['Time - hour of day ',int2str(FFT.jday0)]};
if iihr==2,
	CGRAM.scales{3}=FFT.range(jseq);
	CGRAM.labels{3}='Range - km'
elseif iihr==0,
	CGRAM.scales{3}=FFT.jday0*1440 + FFT.hodax(jseq_out)*60 - FFT.moy0;
	CGRAM.labels{3}=['Time - min'];
end
CGRAM.data=zeros(ntime,nrec,nseq_out);
CGRAM.dlab='Correlation Coefficient';
%
% This is required to get 0 dB at zero time lag for auto correlations:
% Note that nfftx is the length of the IFFT with zero padding or enveloping:`
fft_norm=nfftx/2;
if iienv, fft_norm=2*fft_norm; end
%
ns_out=0;
for js=1:nseq,
	js_out=jseq(js);
	[XH,cspec1]=fft_read1(FFT,js_out,jrec(:,1),jfreq);
	[XH2,cspec2]=fft_read1(FFT2,js_out,jrec(:,2),jfreq);
    if T_exp==-1,
        cspec1=cspec1./max(1.e-20,abs(cspec1)); cspec2=cspec2./max(1.e-20,abs(cspec2));
    elseif T_exp>0,
        disp('not implemented yet');
    end
	if iihann, cspec1=cspec1.*hann_arr; cspec2=cspec2.*hann_arr; end
	cnorm1(:)=sum(zmagsq(cspec1), 1);
	cnorm2(:)=sum(zmagsq(cspec2), 1);
% Norm such that, when rspec = alpha * cspec, max corr. at tau=0 is unity:
	cnorm(:,:)=repmat(sqrt(cnorm1.*cnorm2),[nfreq,1]);
	cspec=cspec1.*conj(cspec2)./cnorm;
% Otherwise, correlations won't be normalized:
	if js==1 & ns_avg>1, cspec_avg=zeros(nfreq,nrec); end
% Average over sequence number (meant for averaging cross_spectra):
	if ns_avg>1,
		cspec_avg(:,:)=cspec_avg(:,:) + cspec;
		if mod(js,ns_avg)==0, cspec(:,:)=cspec_avg/ns_avg; cspec_avg(:,:)=0; end
	end
	if mod(js,ns_avg)==0,
		ns_out=ns_out + 1;
		if ~iienv,
			fft_full(imap1,:)=fft_norm*cspec;
			fft_full(imap2,:)=conj(fft_full(imap1,:));
% Do inverse FFTs:
%%size(pser),size(ifft(fft_full)),keyboard
			pser(:,:)=ifft(fft_full);
			x(:,:)=real(pser(iit,:));
		else
			fft_full(1:nfreq-iih,:)=fft_norm*cspec(iih+1:nfreq,:);
			fft_full(nfftx-iih+1:nfftx,:)=fft_norm*cspec(1:iih,:);
% Do inverse FFTs:
			pser(:,:)=ifft(fft_full);
			x(:,:)=pser(iit,:);
		end
		CGRAM.data(:,:,ns_out)=x(:,:);
%%keyboard
%%		if ~isdata, CGRAM.t_ref(ns_out)=XH.t0; end
	end
%
%%	[t,pser]=myffti(cspec,FFT,tlims,0,fzpad,iienv);
%	keyboard,
    pause(.001);
	if mod(js,10)==1, disp(['Finished sequence ',int2str(js),' of ',int2str(nseq)]); end
end
if iienv==1,
	CGRAM.data=abs(CGRAM.data);
elseif iienv==2,
	CGRAM.data=20*log10(max(1e-30,abs(CGRAM.data)));
elseif iienv==3,
	CGRAM.data=angle(CGRAM.data)/pi;
elseif iienv==4,
	CGRAM.data=unwrap(angle(CGRAM.data))/pi;
end
%
if isempty(clim),
	cmax=max(CGRAM.data(:));
	if iienv==0
		cmin=min(CGRAM.data(:));
	elseif iienv==1,
		cmin=0;
	elseif iienv==2,
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(CGRAM.data(:))/5)*5,cmax-60);
		CGRAM.dlab='Pressure - dB//\muPa^2';
	elseif iienv==3,
		cmin=-1; cmax=1;
		CGRAM.dlab='arg(env)/pi';
	elseif iienv==4,
		cmin=min(CGRAM.data(:));
		CGRAM.dlab='unwrap[arg(env)]/pi';
	elseif iienv==-1,
		cmax=max(dB(CGRAM.data(:)));
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(dB(CGRAM.data(:)))/5)*5,cmax-60);
		CGRAM.dlab='Complex pressure - \muPa';
	end
	clim=[cmin,cmax];
end
%
CGRAM.clim=clim;
CGRAM.fzpad=fzpad; CGRAM.iienv=iienv; 
if iienv==2, CGRAM.iidB=1; elseif iienv==-1, CGRAM.iidB=-1; else, CGRAM.iidB=0; end
%
fclose(FFT.fft_fid);
%
%%keyboard,
if iifig & ~isempty(CGRAM.data),
	view_global; cm=V.cmap; V.cmap=1-gray;
	view5(CGRAM); V.cmap=cm;
end
%
% 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.
