function PRUNE=mcs_prune_init(xyz,f_prune,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=lam_parms(1);
Nlam_xc_max=lam_parms(2);
Nlam_ap_min=lam_parms(3);
Nlam_ap_max=lam_parms(4);
%
[N_ph,d3]=size(xyz);
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);
%
nf_prune=length(f_prune); 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,:));
%
% Compute cumulative apertures in integer units of lambda*Nlam_xc_min:
	cum_ap=round(cum_dist/(lambda(jf)*Nlam_xc_min));
% Divide array into nsub subapertures, consisting of phones
% iisub1(jsub):iisub2(jsub), for jsub=1:nsub;
	iijump=find(diff(cum_ap)>0);
	if length(iijump)<=2, 
		[dum,nq1]=min(abs(cum_dist-.25*aper));
		[dum,nq2]=min(abs(cum_dist-.50*aper));
		[dum,nq3]=min(abs(cum_dist-.75*aper));
		iisub2=unique([nq1;nq2;nq3;nph]);
	else
		iisub2=[iijump ; nph];
	end
	iisub2=[iijump ; nph];
	nsub=length(iisub2);
	iisub1=[1; iisub2(1:nsub-1)+1];
	nph_sub=iisub2 - iisub1 + 1;	% Number of phones in each subaperture
	nph_sub_max=max(nph_sub);		% Max number of phones in subapertures
%
% Make index for placing weighted elements into array bf_mat that can be 
% summed vertically to achieve the subaperture summations:
	idx_bf=zeros(nph,1); sub_cent=zeros(nsub,3);
	for jsub=1:nsub,
		ii=iisub1(jsub):iisub2(jsub);
		sub_cent(jsub,:)=mean(xyz(iiph(ii),:));	% Phase center of subaperture;
		idx_bf(ii)=(jsub-1)*nph_sub_max + [1:nph_sub(jsub)];
	end
%% Example usage to sum weighted elements to get subarray beams:
%%	bf_mat=zeros(nph_sub_max,nsub);
%%	bf_mat(idx_bf)=bf;
%%	bf_sub=sum(bf_mat,1);
%
% Make indices of subapertures to crosscorrelate:
	if nsub>2, nskip=2; else, nskip=1; end
	sub_aper_max=Nlam_xc_max*lambda(jf);	% Max subaperture separation to include;
% Doing it as arrays saves a factor of 10!!:
	lsub=nsub-nskip;
	jsub1=[1:nsub-nskip]' * ones(1,lsub); jsub2=ones(lsub,1)* [1+nskip:nsub];
	sub_comb_all=[1:nsub-nskip]' * [1+nskip:nsub];
	iiok=find(jsub2>=jsub1+nskip);
	delxyz=sub_cent(jsub1(iiok),:) - sub_cent(jsub2(iiok),:);
	sub_aper=sqrt(sum(delxyz.^2,2));
	iiok2=iiok(find(sub_aper<=sub_aper_max));
	idx_xc1=jsub1(iiok2); idx_xc2=jsub2(iiok2); n_xc=length(idx_xc1);
%% Example usage to get correlate pairs of subarray beams and sum 
%% to get PRUNE beamformer output:
%%	bf_mcs=sum(bf_sub(idx_xc1).*bf_sub(idx_xc2));
%
% 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
	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).nph_sub_max=nph_sub_max;
	PRUNE(jf).iisub=[iisub1,iisub2];
	PRUNE(jf).bf_mat=zeros(nph_sub_max,nsub);
	PRUNE(jf).idx_bf=idx_bf;
	PRUNE(jf).n_xc=n_xc;
	PRUNE(jf).idx_xc1=idx_xc1;
	PRUNE(jf).idx_xc2=idx_xc2;
	PRUNE(jf).nsum=nsub*nph_sub_max + n_xc;
%int2str([iisub1,iisub2]),
end
% 
if iifig,
	colors='rbgmcyk'; ncol=length(colors);
	symbols='osdv<>^ph'; nsym=length(symbols);
	for jf=1:nf_prune,
		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);
		for jxc=1:PRUNE(jf).n_xc,
			col=colors(mod(jxc,ncol)+1);
			sym=symbols(mod(jxc,nsym)+1);
			jsub1=PRUNE(jf).idx_xc1(jxc);
			jsub2=PRUNE(jf).idx_xc2(jxc);
%			keyboard
			iisub1=iiph(PRUNE(jf).iisub(jsub1,1):PRUNE(jf).iisub(jsub1,2));
			iisub2=iiph(PRUNE(jf).iisub(jsub2,1):PRUNE(jf).iisub(jsub2,2));
			for jph1=1:length(iisub1),
				plot_cfill(cum_dist0(iisub1(jph1)),cum_dist0(iisub2),'o',col,5);
			end
		end
		axis([0,cum_dist0(N_ph),0,cum_dist0(N_ph)])
		pause(.001);
%		keyboard
	end
end
toc

