function [Am,Bm,Cm,V] = mcmillan(A,B,C,dc);
%
%  mcmillan
%
%  Converts the state-space model A,B,C in damped
%   modal coordinates to Am,Bm,Cm in undamped modal
%   coordinates.
%
%  It is assumed that the states are in complex-conjugate
%  pairs, and that real poles have been removed
%  (e.g. using 'chkpoles')
%
%  'dc' is the vector of eigenvalues of A

%  Reference: Alvin, K.F., "Second-Order Structural 
%   Identification via State Space-Based System
%   Realizations," Ph. D. Dissertation, Dept. of 
%   Aerospace Engineering, Univ. of Colorado, April 1993.

%  Version SWD960205
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

nx = size(A,1);

nm = nx/2;

V = zeros(size(A));

j = sqrt(-1);

%  Define the transformation

for i = 1:nm,

    V([((i-1)*2 + 1):(i*2)],[((i-1)*2 + 1):(i*2)]) = j/(2 * imag(dc(2*i-1))) *...
 		 [dc(2*i),-1;(-1)*dc(2*i-1),1];

end

%  Apply the transformation to the damped-modal model

Vi = inv(V);

Am = real(Vi * A * V);
Bm = real(Vi * B);
Cm = real(C * V);

return
