function PRUNE=mcs_prune_init2(xyz,f_lims,cavg,lam_parms,iifig)
%
% Initializes index vectors and matrices for lambda-subaperture
% matched cross-spectral beamforming.
%
%
% 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.
tic;
Nlam_xc_min=max(.0001,lam_parms(1));
Nlam_xc_max=lam_parms(2);
Nlam_ap_min=max(.0001,lam_parms(3));
Nlam_ap_max=lam_parms(4);
if length(lam_parms)>=5,
	Nph_xc_min=lam_parms(5);
else
	Nph_xc_min=3;
end
%
% Compute pruning frequencies:
dlogf=log10(1.1);		% Allow 10% change between prune frequencies;
logf_lim=log10(f_lims);
nf_prune=max(2,round(diff(logf_lim)/dlogf)),
f_prune=10.^ mylinspace(logf_lim(1),logf_lim(2),nf_prune),
%
[N_ph,d3]=size(xyz); nph=N_ph;
dx=xyz(:,1)*ones(1,nph) - ones(nph,1)*xyz(:,1)';
dy=xyz(:,2)*ones(1,nph) - ones(nph,1)*xyz(:,2)';
dz=xyz(:,3)*ones(1,nph) - ones(nph,1)*xyz(:,3)';
dmat0=sqrt(dx.^2 + dy.^2 + dz.^2);
kph1=zeros(2,100); kph2=kph1; nph_prod=zeros(1,100);
nsub1=zeros(100,1); nsub2=nsub1;
jrow_hi=zeros(50); jcol_lo=jrow_hi;
%
dist0=sqrt(sum(diff(xyz,1).^2,2));
iiuni=all(abs(dist0-dist0(1))./dist0(1)<.01);	% Flag for uniform arrays;
%
cum_dist0=[0; cumsum(dist0)];
aper=cum_dist0(N_ph);
%
lambda=cavg./f_prune;
nsub=zeros(nf_prune,1); 
ii1=1:N_ph-2; ii2=ii1+1;

for jf=1:nf_prune,
% First determine which element indices (iiph) to use at this frequency;
% Discard elements whose spacing is less than Nlam_ap_min apart:
	dist_lam_min=round(cum_dist0/(Nlam_ap_min*lambda(jf)));
	iiok=find(dist_lam_min(2:N_ph-1) > dist_lam_min(1:N_ph-2));
	iiph=[1; 1+iiok; N_ph]; nph=length(iiph);
% If array is longer than Nlam_ap_max wavelengths, prune it:
	if nph<N_ph,
		dist=sqrt(sum(diff(xyz(iiph,:),1).^2,2));
		cum_dist=[0; cumsum(dist)]; 
	else
		dist=dist0; cum_dist=cum_dist0;
	end
	aper=cum_dist(nph); 
	aper_max=Nlam_ap_max*lambda(jf); iileft=0;
	while aper > aper_max,
% Take away phones from end with largest spacing (or alternate if uniform array):
		if iiuni, iileft=~iileft; else, iileft=dist(1)>dist(nph-1); end
		if iileft,
			iiph=iiph(2:nph); dist=dist(2:nph-1); 
			cum_dist=cum_dist(2:nph)-cum_dist(2); nph=nph-1;
		else
			iiph=iiph(1:nph-1); dist=dist(1:nph-2); 
			cum_dist=cum_dist(1:nph-1); nph=nph-1;
		end
		aper=cum_dist(nph);
	end
% Pruning part finished. Compute phase center before doing lambda-apertures:
	phase_cent=mean(xyz(iiph,:));
%
	dmat=dmat0(iiph,iiph);
	xc_min=Nlam_xc_min*lambda(jf); xc_max=Nlam_xc_max*lambda(jf);
%		
	ntri=1; jtri_lo(ntri)=1; jtri_hi(ntri)=nph;
	jtrow0(ntri)=nph; jtcol0(ntri)=1;
	n_xc=0; nph_sub_max=0;
	while ntri>0,
		jtlo=jtri_lo(ntri); jthi=jtri_hi(ntri);
		nrect=1; jcol_lo(nrect)=jtlo; jrow_hi(nrect)=jthi;		
		jtr0=jtrow0(ntri); jtc0=jtcol0(ntri); jtmid=round(.5*(jtr0 + jtc0));
		jrow0(nrect)=jtmid; jcol0(nrect)=jtmid;
% Keep track of min,max rows and cols taken:
		tri_min=[nph nph]; tri_max=[1 1]; n_xc_tri=0;
%%		disp(['Triangle: ',int2str([jtlo jthi jtr0 jtc0 n_xc])]),
	
		while nrect>0,
			jrmax=jrow_hi(nrect); jcmin=jcol_lo(nrect);
			jr0=jrow0(nrect); jc0=jcol0(nrect);
% Find indices of diagonal from idmid,idmid to left and down:
			jrow=[jr0:jrmax]; jcol=[jc0-1:-1:jcmin];
			npos=min(length(jrow),length(jcol)); ii=1:npos;
% Convert 2-D matrix indices to 1-D vector indices:
			j1d=(jcol(ii)-1)*nph + jrow(ii);	
%	
			ap_pos=dmat(j1d);
			iiok=find(ap_pos>=xc_min & ap_pos<=xc_max); nok=length(iiok); 
			nrect=nrect - 1; 
%			
			if nok>0,
				iiok1=iiok(1); iiok2=iiok(nok);
				jrowmax=jrow(iiok2); jcolmin=jcol(iiok2);
				if jrowmax<jrmax,
					jmore=find(dmat(jrowmax:jrmax,jcolmin) <= xc_max);
					jrowmax=jrowmax-1 + jmore(length(jmore));
				end
				if jcolmin>jcmin,
					jmore=find(dmat(jrowmax,jcolmin:-1:jcmin) <= xc_max);
					jcolmin=jcolmin+1 - jmore(length(jmore));
				end
				ns1=jrowmax-jrow(iiok1) + 1; ns2=jcol(iiok1) - jcolmin + 1;
%%		keyboard
				if n_xc==0 | (ns1>=Nph_xc_min & ns2>=Nph_xc_min),
					n_xc=n_xc + 1;
					nsub1(n_xc)=ns1;
% Convert to xyz phone #:
					kph1(1:2,n_xc)=[jrow(iiok1);jrowmax];	
					nsub2(n_xc)=ns2;
					kph2(1:2,n_xc)=[jcolmin;jcol(iiok1)];
					nph_prod(n_xc)=ns1*ns2;
					nph_sub_max=max([nph_sub_max,nsub1(n_xc),nsub2(n_xc)]);
					tri_min=min([tri_min;jrow(iiok1) jcolmin]);
					tri_max=max([tri_max;jrowmax jcol(iiok1)]);
					n_xc_tri=n_xc_tri + 1;
				end
%	ph1=iiph(jrow(iiok1):jrowmax)
%	ph2=iiph(jcolmin:jcol(iiok1))
				if jcolmin>jcmin,
					nrect=nrect + 1;
					jrow0(nrect)=jrow(iiok1); jcol0(nrect)=jcolmin - 1;
					jrow_hi(nrect)=jrowmax; jcol_lo(nrect)=jcmin;
				end
				if jrowmax<jrmax,
					nrect=nrect + 1;
					jrow0(nrect)=jrowmax+1; jcol0(nrect)=jcol(iiok1);
					jrow_hi(nrect)=jrmax; jcol_lo(nrect)=jcolmin;
				end
			end
%%			nrect, keyboard
		end
		ntri=ntri - 1;
		if n_xc_tri>0,
			jtri0=tri_min(1)-1; %jtri0=idmid;
			if jtri0 > jtlo+1,
				ntri=ntri + 1; jtri_lo(ntri)=jtlo; jtri_hi(ntri)=jtri0;
				jtrow0(ntri)=max(1,tri_min(1)-1); 
				jtcol0(ntri)=max(1,tri_min(2)-1);
			end
			jtri0=tri_max(2)+1; %jtri0=idmid;
			if jtri0 < jthi-1,
				ntri=ntri + 1; jtri_lo(ntri)=jtri0; jtri_hi(ntri)=jthi;
				jtrow0(ntri)=min(nph,tri_max(1)+1); 
				jtcol0(ntri)=min(nph,tri_max(2)+1);
			end
		end
	end
% Save quantities in output structure PRUNE:
	if jf==1,
		PRUNE(jf).f_prune=f_prune;	% Include all frequencies in first element
	else
		PRUNE(jf).f_prune=f_prune(jf);
	end
	if n_xc==0,
		PRUNE(1).f_prune(jf)=-999;
	end
	PRUNE(jf).iiph=iiph;
	PRUNE(jf).nph=nph;
	PRUNE(jf).aper=aper;
	PRUNE(jf).lambda=lambda(jf);
	PRUNE(jf).phase_cent=phase_cent;
	PRUNE(jf).cavg=cavg;
	PRUNE(jf).nsub=nsub;
	PRUNE(jf).n_xc=n_xc;
	PRUNE(jf).nph_sub_max=nph_sub_max;
	PRUNE(jf).nsub1=nsub1(1:n_xc);
	PRUNE(jf).nsub2=nsub2(1:n_xc);
	PRUNE(jf).nph_prod=nph_prod(1:n_xc);
	PRUNE(jf).kph1=kph1(1:2,1:n_xc);
	PRUNE(jf).kph2=kph2(1:2,1:n_xc);
	PRUNE(jf).nph_prod_tot=sum(nph_prod(1:n_xc));
% Example usage of PRUNE:
% Compute beamformed weights and sqared magnitudes element-by-element:
%% bf_elements=conj(w) .* X(iiph);
%% bf_norm=real(bf_elements).^2 + imag(bf_elements).^2;
% Form cumulative sums:
%% bf_cum=[0; cumsum( bf_elements )];
%% norm_cum=[0; cumsum( bf_norm )];
% Evaluate subaperture beamformer outputs and norms from cum sums:
%% bf_xc1=bf_cum(kph1(2,:)+1) - bf_cum(kph1(1,:));
%% bf_xc2=bf_cum(kph2(2,:)+1) - bf_cum(kph2(1,:));
%% xc_norm1=norm_cum(kph1(2,:)+1) - norm_cum(kph1(1,:));
%% xc_norm2=norm_cum(kph2(2,:)+1) - norm_cum(kph2(1,:));
% Compute total norm for correlation:
%% xc_norm_tot=sqrt(xc_norm1 .* xc_norm2);
% Crosscorrelate sub-ap BF outputs and sum up (include weighting according to # phones):
%% mcs_bf=sum((bf_xc1 .* bf_xc2) .* nph_prod ./ xc_norm_tot ) / nph_prod_tot;
%% 
%% 
%%	PRUNE(jf).bf_mat=zeros(nph_sub_max,nsub);
	if iifig,
		colors='rbgmcyk'; ncol=length(colors);
		symbols='osdv<>^ph'; nsym=length(symbols);
		figure;
		iiph=PRUNE(jf).iiph;
		plot_cfill(cum_dist0(iiph),cum_dist0(iiph),'s','k',8); hold on; axis('ij');
		col=colors(mod(jf,ncol)+1);
		nsub1=PRUNE(jf).nsub1; nsub2=PRUNE(jf).nsub2;
		for jxc=1:PRUNE(jf).n_xc,
			col=colors(mod(jxc,ncol)+1);
			sym=symbols(mod(jxc,nsym)+1);
			kph1x=PRUNE(jf).kph1(1:2,jxc); kph2x=PRUNE(jf).kph2(1:2,jxc);
			iiph1=iiph(kph1x(1):kph1x(2)); iiph2=iiph(kph2x(1):kph2x(2));
			for jph1=1:length(iiph1),
				plot_cfill(cum_dist0(iiph2),cum_dist0(iiph1(jph1)),'o',col,5);
			end
		end
		axis([0,cum_dist0(N_ph),0,cum_dist0(N_ph)])
		pause(.001); %%keyboard
	end
end
% 
toc
