function [BEAMS,REP]=mcs_rep_init(beam_parm)
%
% Computes replica vectors for array geometry given in xyz at frequncies
% in f_rep.  Outputs beamforming information in structure BEAMS, replica vectors
% in REP.
%
% beam_parm	= [iisin,iivar,fac_var], where iisin=1 for sin space, iivar=1 for
%			  nphi_bf dependent on frequency, fac_var is beam spacing factor
%			  used when iivar=1 [nphi_bf=fac_var*aper/(lambda/2)];
%
%
% 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.
mcs_global
%
% Compute replica vector frequencies for frequencies specified in MCS.freq_lim:
nbinh=floor(MCS.nbin/2); nbinh2=MCS.nbin - nbinh - 1; df=FFT.df;
[nr,nc]=size(MCS.freq_lim),
if nc==1,
	MCS.f_rep=MCS.freq_lim;
elseif nr==1 & nc==2,
	df_rep=MCS.nbin*df;
% Make replica vector center frequencies coincide with bin frequencies:
	MCS.f_rep=[MCS.freq_lim(1)+nbinh*df:df_rep:MCS.freq_lim(2)-nbinh2*df]';
else
	error('freq_lim bad');
end
%
xyz=MCS.xyz; cavg=MCS.cavg; f_rep=MCS.f_rep; phi_lim=MCS.phi_lim; th_lim=MCS.th_lim;
% Read off options associated with variable numbers of beams vs. frequency:
tic
iisin=beam_parm(1); iivar=beam_parm(2); 
if iivar, fac_var=beam_parm(3); end
% Read off pruning frequencies from first element of structure PRUNE:
f_prune=PRUNE(1).f_prune;
[N_ph,dum]=size(xyz);
%
nf_rep=length(f_rep); lambda=cavg./f_rep;
% Pre-compute xyz positions of all elements in the array shifted to phase center:
phase_cent=PRUNE(1).phase_cent;		% Use a single phase center for all freqs:
xyz_cent=xyz - ones(N_ph,1)*[phase_cent(1:2),0];	% Don't change z values;
% Compute LOB in deg E of N:
pifac=180/pi;
iiph=PRUNE(1).iiph; nph=PRUNE(1).nph;
pfit=polyfit(xyz(iiph,1),xyz(iiph,2),1);
if pfit(1)~=0, lob1=atan(1/pfit(1)); else, lob1=pi/2; end
delxy=diff(xyz([iiph(1),iiph(nph)],1:2),1,1);
lob2=atan2(delxy(1),delxy(2));
if abs(lob2-lob1)>pi/2, 
	if lob1>0, lob=lob1-pi; else, lob=lob1+pi; end
else
	lob=lob1;
end
lob_deg=lob*pifac;
if iisin,
% For beams uniformly spaced in sin space relative to LOB, rotate xyz so LOB=90 deg:
	rot_mat=[sin(lob) -cos(lob);cos(lob) sin(lob)];
	xyz_rot=[xyz_cent(:,1:2)*rot_mat , xyz_cent(:,3)];
%%	delxy=diff(xyz_rot([iiph(1),iiph(nph)],1:2),1,1);
%%	lob_check=atan2(delxy(1),delxy(2))*180/pi;
else
	xyz_rot=xyz_cent;
end
%keyboard,
%
% Pre-compute direction cosines associated with steering vectors:
pirad=pi/180;
theta=adi_lim_eval(th_lim);
ntheta=length(theta); costheta=cos(theta(:)'*pirad); sintheta=sin(theta(:)'*pirad);
if ~iivar,
	phi_bf=adi_lim_eval(phi_lim)'; nphi_bf=length(phi_bf); nbeam=ntheta*nphi_bf;
	phi_out=phi_bf; nphi_out=nphi_bf;
%%	if iisin, phi_bf=(phi_bf + 1)*90 - 90; end % NO, this give uniform in phi!!
	[sinphi,cosphi]=sin_cos_calc(phi_out,iisin);
	dcx=sinphi*costheta; dcy=cosphi*costheta; dcz=ones(nphi_bf,1)*sintheta;
	dc_rep=[dcx(:)';dcy(:)';dcz(:)'];	% 3 by nbeam matrix of direction cosines
else
% For nphi_bf frequency-dependent, wait to compute nphi_bf, phi_bf until inside freq loop:
	phi_out=adi_lim_eval(phi_lim); nphi_out=length(phi_out);
% Make sure phi_out is uniformly spaced:
	phi_out=mylinspace(phi_out(1),phi_out(nphi_out),nphi_out);
	if iisin, [sinphi,cosphi]=sin_cos_calc(phi_out,iisin); end
	idx_st=zeros(3,nphi_out); qfac_st=zeros(3,nphi_out);
end
if iisin,
	phi_sin=phi_out; phi_deon=atan2(sinphi,cosphi)*180/pi + (lob_deg-90);
else
	phi_deon=phi_out; phi_sin=phi_out;
% phi_sin is not sin space, but for this case we do not care about it!
end
% Compute wavenumber xk=2*pi*f/c:
xk_rep=(2*pi/cavg)*f_rep;
% Positive for replica vector (ABF), negative for steering vector (CBF)
%
if MCS.abf_parm(1)>0, iiabf=1; sg=1; else, iiabf=0; sg=-1; end
%
freq_abf=zeros(1,nf_rep)
nphmax=0; nbeammax=0;
jf_prune=zeros(nf_rep,1);
jbin=zeros(nf_rep,2);
% Loop over replica vector frequencies:
for jf=1:nf_rep,
	freq=f_rep(jf); xk=xk_rep(jf);
% Find closest pruning frequency and find phone numbers used from PRUNE:
	[df,jfp]=min(abs(f_rep(jf)-f_prune));
	jf_prune(jf)=jfp; aper=PRUNE(jfp).aper; 
	iiph=PRUNE(jfp).iiph; nph=PRUNE(jfp).nph; aper=PRUNE(jfp).aper;
	nphmax=max(nphmax,nph);
	xk_xyz=xk*xyz_rot(iiph,:);
	if iivar,
		hlam=.5*PRUNE(jfp).lambda;
		nphi_bf=max(2,round(fac_var*aper/hlam)); nbeam=ntheta*nphi_bf;
		phi_bf=mylinspace(phi_lim(1),phi_lim(2),nphi_bf);
		[sinphi,cosphi]=sin_cos_calc(phi_bf,iisin);
		dcx=sinphi*costheta; dcy=cosphi*costheta; dcz=ones(nphi_bf,1)*sintheta;
		dc_rep=[dcx(:)';dcy(:)';dcz(:)'];	% 3 by nbeam matrix of direction cosines
% Compute quadratic factors for interpolating the nphi_bf steered beams in phi_bf 
% to the nphi_out output beams:
% Original quadratic scheme left discontinuities across phi_bf midpoints because
% two different quadratic fits gave different values.  So average the two fits:
		phi_out=mylinspace(phi_lim(1),phi_lim(2),nphi_out);
		dphi_bf=(phi_bf(nphi_bf) - phi_bf(1))/max(1,nphi_bf-1);
%		idx_st(1,:)=round((phi_out - phi_bf(1))/dphi_bf) + 1;
%		idx_st(2,:)=max(idx_st(1,:) - 1,1);
%		idx_st(3,:)=min(idx_st(1,:) + 1,nphi_bf);
%		x=(phi_out - phi_bf(idx_st(1,:)))/dphi_bf; xsq=x.^2;
%		qfac_st(1,:)=1 - xsq;		% 3 by nphi_out
%		qfac_st(2,:)=.5*(xsq - x);
%		qfac_st(3,:)=.5*(xsq + x);
		idx_st(2,:)=floor((phi_out - phi_bf(1))/dphi_bf) + 1;
		idx_st(1,:)=max(idx_st(2,:) - 1,1);
		idx_st(3,:)=min(idx_st(2,:) + 1,nphi_bf);
		idx_st(4,:)=min(idx_st(2,:) + 2,nphi_bf);
		x1=(phi_out - phi_bf(idx_st(2,:)))/dphi_bf; x1sq=x1.^2;
		x2=(phi_out - phi_bf(idx_st(3,:)))/dphi_bf; x2sq=x2.^2;
		qfac_st(1,:)=.5 * (.5*(x1sq - x1) );
		qfac_st(2,:)=.5 * (1 - x1sq + .5*(x2sq - x2) );		% 4 by nphi_out
		qfac_st(3,:)=.5 * (.5*(x1sq + x1) + 1 - x2sq );
		qfac_st(4,:)=.5 * (.5*(x2sq + x2));
% Example usage of idx_st and qfac_st for interpolation from bf_st to bf_out:
%% bf_out=sum( bf_st(idx_st) .* qfac_st , 1);
	end
	nbeammax=max(nbeammax,nbeam);
	phase=xk_xyz(:,1)*dc_rep(1,:) + xk_xyz(:,2)*dc_rep(2,:) + xk_xyz(:,3)*dc_rep(3,:);
	REP(jf).type='PW';
	REP(jf).beam_parm=beam_parm;
	REP(jf).f_rep=f_rep(jf);
	REP(jf).jf_prune=jfp;
	REP(jf).nbeam=nbeam; REP(jf).phi_bf=phi_bf; REP(jf).ntheta=theta;
	REP(jf).phi_bf=phi_bf;
	REP(jf).theta=theta;
%
	if iiabf,
		REP(jf).repvec=cos(phase) + i*sin(phase);	% nph by nbeam matrix;
		REP(jf).wtsvec=zeros(nph,nbeam);			% nph by nbeam matrix;
		REP(jf).R=eye(nph,nph) + (eps*i)*eye(nph,nph);	% nph by nph CSDM;
		freq_abf(jf)=1;
% Gw_ref is the reference white noise gain: 10*log10(N) for f-ind WNGC, 
% 10*log10(min(N,L/(lambda/2))) for f-dep WNGC.
		if MCS.abf_parm(5),
			REP(jf).Gw_ref=10*log10(min(nph,aper/(lambda(jf)/2))); 
		else
			REP(jf).Gw_ref=10*log10(nph); 
		end
	else
		REP(jf).wtsvec=(1/nph)*(cos(phase) - i*sin(phase));	% nph by nbeam matrix
	end
	if iivar,
		REP(jf).idx_st=idx_st; REP(jf).qfac_st=qfac_st;
	end
end
%
BEAMS.type='PW'; BEAMS.iivar=iivar;
BEAMS.beam_parm=beam_parm;
BEAMS.nbeammax=nbeammax;
BEAMS.nf_rep=nf_rep;
BEAMS.f_rep=f_rep;
BEAMS.df_rep=(f_rep(nf_rep) - f_rep(1))/max(1,nf_rep-1);
BEAMS.jf_prune=jf_prune;
BEAMS.lob=lob; BEAMS.lob_deg=lob_deg;
BEAMS.nbeam=nphi_out*ntheta; BEAMS.nphi_out=nphi_out; BEAMS.ntheta=ntheta;
BEAMS.phi_out=phi_out; BEAMS.phi_deon=phi_deon; BEAMS.phi_sin=phi_sin;
BEAMS.nphmax=nphmax;
BEAMS.jf_abf=find(freq_abf); BEAMS.nf_abf=length(BEAMS.jf_abf);
%%%
BEAMS.nf_rep=length(MCS.f_rep);
% Find FFT bins associated with each replica vector center frequency:
df=FFT.df; fax=(FFT.imap(:,1) - 1) * df;
nbinh=floor(MCS.nbin/2); nbinh2=MCS.nbin - nbinh - 1;
iif_cent=round((MCS.f_rep-fax(1))/df) + 1;
BEAMS.jf_rep2bin=[iif_cent - nbinh,iif_cent + nbinh2];	% nf_rep by 2
% Find FFT bins that need to be read:
%%BEAMS.iifbin=[BEAMS.jf_rep2bin(1,1) : BEAMS.jf_rep2bin(BEAMS.nf_rep,2)];
% Make life easier by reading in frequencies up to last one needed:
BEAMS.iifbin=[1 : BEAMS.jf_rep2bin(BEAMS.nf_rep,2)];
% Find frequencies that will be beamformed:
nf_bf=0; jbin2p=0; nbinmax=0;
BEAMS.f_bf=zeros(2*MCS.nbin*BEAMS.nf_rep,1); BEAMS.jf_bin2rep=BEAMS.f_bf;
%%BEAMS.iifbin=BEAMS.f_bf;
for jf_rep=1:BEAMS.nf_rep,
	jbin1=max(jbin2p+1,BEAMS.jf_rep2bin(jf_rep,1)); 
	jbin2=BEAMS.jf_rep2bin(jf_rep,2);
	nbf=jbin2-jbin1+1; iif=nf_bf+1:nf_bf + nbf; nbinmax=max(nbinmax,nbf);
	BEAMS.f_bf(iif)=fax(jbin1:jbin2);
	BEAMS.jf_bin2rep(iif)=jf_rep;
%%	BEAMS.iifbin(iif)=[jbin1:jbin2];
	jbin2p=jbin2; nf_bf=nf_bf + nbf;
end
BEAMS.nbinmax=nbinmax;
BEAMS.nf_bf=nf_bf;
BEAMS.f_bf=BEAMS.f_bf(1:nf_bf); 
BEAMS.jf_bin2rep=BEAMS.jf_bin2rep(1:nf_bf); 
%%BEAMS.iifbin=BEAMS.iifbin(1:nf_bf);
%
%
% Now map bin frequencies in FFT file to replica vectors:
% Find number of FFT bins to be associated with each replica vector:
%nbinh=floor(MCS.nbin/2); nbinh2=nbinh+MCS.nbin-1;
%fax=(FFT.imap(:,1) - 1) * FFT.df;
%[nr,nc]=size(MCS.freq_lim);
%if nc==1,
% For a list of center frequencies that should not overlap:
%	iif_cent=round((MCS.f_rep-fax(1))/FFT.df) + 1;
%	nf_bf=0; BEAMS.iifbin=zeros(1,MCS.nbin*BEAMS.nf_rep);
%	for jf=1:BEAMS.nf_rep,
%		jbin1=max(1,iif_cent(jf)-nbinh); 
%		jbin2=min(FFT.nfreq,iif_cent(jf)-nbinh2);
%		nf_bf0=jbin2-jbin1 + 1;
%		BEAMS.iifbin(nf_bf+1:nf_bf + nf_bf0)=[jbin1:jbin2];
%		nf_bf=nf_bf + nf_bf0;
%	end
%	BEAMS.iifbin=BEAMS.iifbin(1:nf_bf);
%elseif nr==1 & nc==2,
% For a continuous frequency band:
%	BEAMS.iifbin=find(fax>=MCS.freq_lim(1)-FFT.df/2 & fax<=MCS.freq_lim(2)+FFT.df/2);
%else
%	error('freq_lim bad');
%end
%BEAMS.f_bf=fax(BEAMS.iifbin);
%BEAMS.nf_bf=length(BEAMS.iifbin);

toc
%%%
%
function [sinphi,cosphi]=sin_cos_calc(phi_bf,iisin);
%
if iisin,
	sinphi=phi_bf(:);
	ii1=find(sinphi>1 & sinphi<=3); 	sinphi(ii1)=2-sinphi(ii1); 
	ii3=find(sinphi>3); 				sinphi(ii3)=sinphi(ii3)-4;
	iim1=find(sinphi<-1 & sinphi>=-3); 	sinphi(iim1)=-2-sinphi(iim1);
	iim3=find(sinphi<-3); 				sinphi(iim3)=sinphi(iim3)+4;
	cosphi=sqrt(1 - sinphi.^2);
	cosphi(ii1)=-cosphi(ii1);
	cosphi(iim1)=-cosphi(iim1);
else
	cosphi=cos(phi_bf(:)*pirad); sinphi=sin(phi_bf(:)*pirad);
end
