function [A,B,C,D] = modalLTI(wn,z,m,r,randomized,prd,dt)
% [A,B,C,D] = modalLTI(wn,z,m,r,randomized,prd.dt)
% Generate a state-space realization for a random 
% strictly-proper multi-input multi-output
% linear time-invariant system in its diagonalized (modal) basis. 
% The eigenvectors of the diagonalized system are linearly independent, 
% so repeated eigenvalues are allowed.   
%
%    INPUT              DESCRIPTION
%    =====              ============================
%     wn                natural frequencies of each mode  (1xn) , rad/s
%     z                 damping ratio for each mode (0<z<0.99) default:  0.10  
%     m                 number of outputs
%     r                 number of inputs
%     randomized        choose 'rand' for a random B and C    default: 'rand'
%     prd               choose 'prd' for positive-real dissipative system  
%     dt                time step for discrete-time systems (optional)
%
%    OUTPUT             DESCRIPTION
%    =====              ============================
%    A,B,C,D            a diagonalized state-space realization 
%
% Examples
% ========
%  [A.B.C.D] = modalLTI( 2*pi*[1:10], 0.05, 4, 3); 
%  [A.B.C.D] = modalLTI( 2*pi*10*rand(1,10),0.05*rand(1,10),4,3,'rand','prd'); 
%
%   System characteristics spcified by algorithm choices:
%
%          |  randomized B, C         |  not randomized B, C
% ---------+--------------------------+---------------------
%     prd  |  positive real           |  positive real
% not prd  |  non-minimum-phase zeros |  minimum-phase zeros
%
%  H.P. Gavin, Duke Univ., 2017-09-28

 if nargin < 1, wn = 2*pi*linspace(1,30,5);  end;    
 if nargin < 2, z  = 0.10*ones(size(wn)); end;    
 if nargin < 3, m  = 3; end
 if nargin < 4, r  = 3; end
 if nargin < 5, randomized = 'xunif'; end
 if nargin < 6, prd = 'xprd'; end

 z = min(max(z,0),0.99);

 rr = max(m,r); mm = rr; % Pos.Real and Passive  systems are square

 n = 2*length(wn);

 i = sqrt(-1.0);

 wn = sort(wn);
 wn = wn(:); z = z(:);
 eigA =  wn .* (sqrt(z.^2-1)-z);

 switch randomized  

   case 'rand'
     b  = rand(n/2,rr) + i*rand(n/2,rr);
     c  = rand(mm,n/2) + i*rand(mm,n/2);

   case 'unif'
     idx = [1:n/2]'*[1:rr];  
     b  = idx + i*(idx);
     c  = 1   + i*(idx');

   otherwise                  %  Positive Real
     q = diag([1:n/2])*ones(n/2,1)*[1:rr].^2*pi/n;
     p = [1:mm]'.^2*ones(1,n/2)*diag([1:n/2])*pi/n;
     b =   cos( q ) + i*sin( q );
     c  =  cos( p ) - i*sin( p );

 end

 pol = zeros(1,n);
 pol( 1:2:n )  = eigA; pol( 2:2:n) = conj(eigA);   A = diag(pol);
 B( 1:2:n , : ) = b;   B( 2:2:n , :  ) = conj(b);
 C( : , 1:2:n ) = c;   C(  : , 2:2:n ) = conj(c);
 D  = zeros(m,r);                       % strictly proper

%cb_err_err = (C-B') ./ C

 if strcmp(prd,'prd'),
   B = abs(real(B));  % make the system dissipative
   C = abs(real(C));
   C = B';            % simple positive-real lemma for D=0 and P=I
 end    
 B = B(:,1:r);
 C = C(1:m,:);

 if nargout == 0 % testing ...
   damp(A)
   mimoBode(A,B,C,D,2*pi*logspace(-0.5,2,500),[],93,'b');
%  P = liap(A,0.001*ones(n,rr));
%  nothin = find(P-P')
%  neg_def = eig(P*A+A'*P)
%  pos_def = eig(P)
   dT = 0.005;
   Bd = A \ ( expm(A*dT) - eye(n) ) * B;
   Ad = expm(A*dT); 
   if ( n ~= rank(ctrb(Ad,Bd,50)) ), disp(' (A,B) not controllable '); end
   if ( n ~= rank(obsv(Ad,C ,50)) ), disp(' (A,C) not observable '); end
   nT = 400; dT = 0.005; t = dT*[0:nT-1]; u = 1*randn(r,nT)/sqrt(dT);
   y = lsym(A,B,C,D,u,t,1*ones(n,1));
   figure(95)
   plot(t,y);
   power = (u*y')*dT/nT

   zro = invzeros(A,B,C,D) 
   pz_plot(pol,zro);

 end

 if nargin == 7          % ... continuous to discrete time transformation
   B = A \ ( expm(A*dt) - eye(n) ) * B;
   A = expm(A*dt); 
 end

% =======================================================================
% revised: 2017-09-28, 2017-12-11

