function mpc = mpc_comp(C);
%
%  MPC_COMP
%
%  mpc = mpc_comp(C);
%
%  Computes the modal phase collinearity (MPC) for
%  the complex mode shapes in C (damped modal coordinates)
%  It is assumed that C consists of complex-conjugate pairs
%  of modes (i.e. real poles have been removed)

%  Reference: "Consistent-Mode Indicator for the Eigensystem
%    Realization Algorithm," Pappa, Elliott, Schenk, J. of
%    Guidance, Control, and Dynamics, v. 16, No. 5, 1993.

%  Version SWD960206
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This matlab source code was originally     %
% developed as part of "DIAMOND" at          %
% Los Alamos National Laboratory. It may     %
% be copied, modified, and distributed in    %
% any form, provided:                        %
%  a) This notice accompanies the files and  %
%     appears near the top of all source     %
%     code files.                            %
%  b) No payment or commercial services are  %
%     received in exchange for the code.     %
%                                            %
% Original copyright is reserved by the      %
% Regents of the University of California,   %
% in addition to Scott W. Doebling, Phillip  %
% J. Cornwell, Erik G. Straser, and Charles  %
% R. Farrar.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mpc_out = zeros(size(C,2),1);

for i = 1:2:size(C,2);

  mode = C(:,i);
  rm = real(mode);
  im = imag(mode);

  Sxx = rm' * rm;
  Syy = im' * im;
  Sxy = rm' * im;

  eta = (Syy - Sxx) / (2 * Sxy);

  lam1 = (Sxx + Syy)/2 + Sxy * sqrt(eta^2 + 1);
  lam2 = (Sxx + Syy)/2 - Sxy * sqrt(eta^2 + 1);

  mpc(i:i+1) = ones(2,1) * ((lam1-lam2)/(lam1 + lam2))^2;

end

mpc = mpc';

return
