function mcs_csdm_upd
%
%
% 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
%
alpha=MCS.alpha; alph1=MCS.alph1;
R_cur=zeros(BEAMS.nphmax,BEAMS.nphmax);
U=R_cur; D=R_cur;
%
%%X=ones(BEAMS.nbinmax,BEAMS.nphmax);
% Form current CSDM and update time-averaged CSDM:
for jf=1:BEAMS.nf_abf,
	jf_rep=BEAMS.jf_abf(jf);
	iibin=BEAMS.jf_rep2bin(jf_ref,1):BEAMS.jf_rep2bin(jf_ref,2);
	jfp=REP(jf_rep).jf_prune; ii=1:PRUNE(jfp).nph;
	iiph=PRUNE(jfp).iiph;
%%	X(iif,ii)=xspec(iibin,iiph);	% nfbin_avg by nph;
	R_cur(ii,ii)=xspec(iibin,iiph)' * xspec(iibin,iiph);	% nph by nph (avg over nfbin_avg);
	REP(jf_rep).R = alph1*REP(jf_rep).R + alpha*R_cur(ii,ii);
	[U(ii,ii),D(ii,ii)]=eig(R);
%
	nbeam=REP(jf_rep).nbeam;
	evals=diag(D)'; emax=max(evals); evals=evals/emax;
	if ii_mvdr,
% Full rank MVDR (repvec is nph by nbeam):
		evals_mat=1./(evals.'*ones(1,nbeam));
		Uhv_mat=U'*REP(jf).repvec; % (nph,nbeam)
		vuuv_mat=real(Uhv_mat).^2 + imag(Uhv_mat).^2; % (nph,nbeam)
		den_mat=ones(nph,1)*sum(evals_mat.*vuuv_mat);
%
		Uhv_diag=evals_mat.*Uhv_mat;
		num_wts=U*Uhv_diag;
		REP(jf_rep).wtsvec=num_wts./den_mat; % (nph,nbeam)
%
%%		P_out=real(sum(conj(wts_dmr).*(R*wts_dmr)));
	else
% DMR:
		eps_dmr=MCS.abf_parm(3);
		D=min(nph,MCS.abf_parm(2)); iiD=1:D;
		nph_vec=nph*ones(1,nbeam);
% Set up (nph,nbeam) matrices:
		evals_mat=evals(iiD).'*ones(1,nbeam);
		eps_mat=ones(D,1)*eps_dmr;
		beta_mat=evals_mat./(evals_mat + eps_mat); % (D,nbeam)
% Compute denominator of weights (repvec is nph by nbeam):
		Uhv_mat=U(:,iiD)'*REP(jf).repvec; % (D,nbeam)
		vuuv_mat=abs(Uhv_mat).^2; % (D,nbeam)
		den_summand=beta_mat.*vuuv_mat;
		den=nph_vec - sum(den_summand); %(1,nbeam)
		den_mat=ones(nph,1)*den;  %(nph,nbeam)
%
		Uhv_diag=beta_mat.*Uhv_mat;
		num_wts=REP(jf).repvec - U(:,iiD)*Uhv_diag;
		wts_dmr=num_wts./den_mat; % (nph,nbeam)
%
		P_out=real(sum(conj(wts_dmr).*(R*wts_dmr)));
	end

end

	
