function zz = syszero(A,B,C,D,tol)
% zeros = syszero(A,B,C,D,tol)
% invariant and decoupling zeros of a continuous-time LTI system
% the condition for an invariant zero is that the pencil [ zI-A, -B; C D ]  is
% rank deficient.  For zeros that are not poles (i.e., for minimal realizaitons)
% invariant zeros, z, make H(z) rank-deficient.
% method: use QZ decomposition to find the general eigenvalues of the
% Rosenbrock system matrix,  padded with zero or randn to make it square.
% tol : tolerance value for rank determination,  default = 1e-4

 [nn,rr,mm] = abcddim(A,B,C,D);
  
 if nargin < 5, tol = 1e-6; end

% make the system square by padding with zeros or randn, as needed (cluge?)
  re = mm-rr;   % rank(C) - rank(B);
  me = rr-mm;   % rank(B) - rank(C);
  rm = max(mm,rr);

  zivi = [];        % zeros and associated eigenvectors
  for iter = 1:4

    Bx = B;  Cx = C;  Dx = D;
    if iter == 1  % zero padding for system zeros 
      if mm > rr, Bx = [ B , zeros(nn,re) ]; Dx = [ D , zeros(mm,re) ]; end
      if mm < rr, Cx = [ C ; zeros(me,nn) ]; Dx = [ D ; zeros(me,rr) ]; end
    else          % randn padding for invariant zeros
      if mm > rr, Bx = [ B , randn(nn,re) ]; Dx = [ D , randn(mm,re) ]; end
      if mm < rr, Cx = [ C ; randn(me,nn) ]; Dx = [ D ; randn(me,rr) ]; end
    end

    abcd = [ -A , -Bx ; Cx , Dx ];	% Rosenbrock System Matrix
    I0 = [ eye(nn) , zeros(nn,rm) ; zeros(rm,nn+rm) ];
    [vv,zz] = eig ( abcd , I0, 'qz' );
    zz = -diag(zz);
    idx = isfinite(zz);
    if iter == 1
      z1v1 = [ zz(idx)' ; vv(1:nn+rr,idx) ];
    else
      zivi = [ zivi , [ zz(idx)' ; vv(1:nn+rr,idx) ] ];
    end
  end
 
 zv = unique_cols( [z1v1 zivi] , tol );

 zz = zv(1,:);

 nz  = length(zz);

% Are both the pencil of the Rosenbrock System Matrix _and_
% the transfer function matrix rank deficient??
% confirm that all the zeros are invariant zeros

 nrcABCD = min(size([ A B ; C D ]));
 min_mr = min(size(D));
 pp = eig(A);

 good_zero_index = [];
 fprintf('                   rank  rank  rank  rank \n')
 fprintf('     zero         zABCD     G   zAB   zAC sys inv trans inp outp\n')
 fprintf('  ______________________________________________________________\n')
 for jj = 1:nz
   rank_G = min_mr;
   rank_zAB   = rank( [ zz(jj)*eye(nn) - A , -B ] , tol );  %  inp decpl zero
   rank_zAC   = rank( [ zz(jj)*eye(nn) - A ;  C ] , tol );  % outp decpl zero
   rank_zABCD = rank( [ zz(jj)*eye(nn) - A , -B ; C , D ] , tol ); % inv zero
%  cn_zABCD   = cond( [ zz(jj)*eye(nn) - A , -B ; C , D ] );   % sanity check
   if all(pp ~= zz(jj))
     rank_G = rank(C*inv(zz(jj)*eye(nn)-A)*B + D , tol );
   end
   if (rank_zAB<nn) || (rank_zAC<nn) || (rank_zABCD<nrcABCD) || (rank_G<min_mr)
     good_zero_index = [ good_zero_index , jj ];

     fprintf(' %15.6f %6d %5d %5d %5d', zz(jj), rank_zABCD, rank_G, rank_zAB, rank_zAC );

     fprintf('   *');

     if ( rank_zABCD < nrcABCD) fprintf('   *'); else fprintf('    ');   end

     if ( rank_G < min_mr ) fprintf('     *');   else fprintf('      '); end

     if ( rank_zAB < nn ) fprintf('   *');       else fprintf('    ');   end

     if ( rank_zAC < nn ) fprintf('    *');      else fprintf('     ');  end

     fprintf('\n');
   end
       
 end
 zz  =  zz(good_zero_index);       % return only the invariant zeros

end
%---------------------------------------------------------------------
%   modified from the .m-function tzero2.m by A.S. Hodel, August 1993 
%
%   References 
%
%   Hodel, Computation of Zeros with Balancing, 1992 Lin. Alg. Appl.
%
%   Emami-Naeini and Van Dooren, Automatica, 1982.
%
%   A.J.G. MacFarlane and N. Karcanias,
%   Poles and zeros of linear multivariable systems:
%    A survey of the algebraic, geometric and complex-variable theory
%   International Journal of Control, 24(1):33-74, July 1976
%
%   from : R. Bruce Tenison <btenison@eng.auburn.edu>
%   Created: July 4, 1994
%   A. S. Hodel Aug 1995: allow for MIMO and system data structures


function X = unique_cols( X, tol )
% X = unique_cols( X, tol )
% eliminate repeated columns of X
% Columns are classified as being equal to each other
% within a tolerance, tol, if the absolute value of their differences 
% is less than tol, for all elements of the columns

 [m,n] = size(X);

 idx = ones(1,n);
 for jj = 1:n
   for kk = jj+1:n 
     if all( abs ( X(:,jj) - X(:,kk) ) < tol ) 
       idx(kk) = 0;
     end
   end
 end

 X = X(:,find(idx));

end
% H.P. Gavin,  2017-11-09 , 2022-12-28
% ----------------------------------------------------------------------------
