% function out = szeros(sys,epp)
%
%   Finds the transmission zeros of a SYSTEM matrix. 
%   Occasionally, large zeros are included which may 
%   actually be at infinity. Solving for the transmission 
%   zeros of a system involves two generalized eigenvalue 
%   problems. EPP (optional) defines if the difference between 
%   two generalize eigenvalues is small.
%
%   See also: EIG and SPOLES. 

%   Algorithm based on Laub & Moore 1978 paper, Automatica

function out = szeros(sys,epp)
  if nargin < 1
   disp('usage: szeros(sys)')
   return
  end
  if nargin == 1
   epp = eps;
  end
  [mtype,mrows,mcols,mnum] = minfo(sys);
  if mtype == 'syst'
    [a,b,c,d] = unpck(sys);
    if mnum == 0
      disp('SYSTEM has no states')
    end
    sysu = [a b; c d];

% find generalized eigenvalues of a square system matrix

    if mrows == mcols
      x = zeros(mnum+mcols,mnum+mcols);
      x(1:mnum,1:mnum) = eye(mnum);
      z = eig(sysu,x);
      out = z(~isnan(z) & finite(z));
    else
      nrm = norm(sysu,1);
      if mrows > mcols
        x1 = [ sysu  nrm*(rand(mnum+mrows,mrows-mcols)-.5)];
        x2 = [ sysu  nrm*(rand(mnum+mrows,mrows-mcols)-.5)];
      else
        x1 = [ sysu; nrm*(rand(mcols-mrows,mnum+mcols)-.5)];
        x2 = [ sysu; nrm*(rand(mcols-mrows,mnum+mcols)-.5)];
      end

      x = zeros(x1);
      x(1:mnum,1:mnum) = eye(mnum);
      z1 = eig(x1,x);
      z2 = eig(x2,x);
      z1 = z1(~isnan(z1) & finite(z1));
      z2 = z2(~isnan(z2) & finite(z2));
      nz = length(z1);
      for i=1:nz
        if any(abs(z1(i)-z2) < nrm*sqrt(epp))
          out = [out; z1(i)];
        end
      end
    end
  else
     error('matrix is not a SYSTEM matrix')
     return
  end 
%
% Copyright MUSYN INC 1991,  All Rights Reserved
