function [z,p,k] = ss2zp(a,b,c,d,iu)
%SS2ZP  State-space to zero-pole conversion.
%   [Z,P,K] = SS2ZP(A,B,C,D,IU)  calculates the transfer function in
%   factored form:
%
%                     -1          (s-z1)(s-z2)...(s-zn)
%       H(s) = C(sI-A) B + D =  k ---------------------
%                                 (s-p1)(s-p2)...(s-pn)
%   of the system:
%       .
%       x = Ax + Bu
%       y = Cx + Du
%
%   from the single input IU.  The vector P contains the pole 
%   locations of the denominator of the transfer function.  The 
%   numerator zeros are returned in the columns of matrix Z with as 
%   many columns as there are outputs y.  The gains for each numerator
%   transfer function are returned in column vector K.
%
%   See also ZP2SS,PZMAP,TZERO, EIG.

%   Copyright 1984-2014 The MathWorks, Inc.
narginchk(4,5)
[msg,a,b,c,d]=abcdchk(a,b,c,d); error(msg);
[nx,~] = size(a);

if nargin<5
   if nx>0,
      [~,nu] = size(b);
   else
      [~,nu] = size(d);
   end
   if (nu<=1),
      iu = 1;
   else
      error(message('Controllib:general:NeedIU'))
   end
end

% Remove relevant input:
if ~isempty(b), b = b(:,iu); end
if ~isempty(d), d = d(:,iu); end

% Do poles first
p = eig(a);

% Compute zeros and gains
[ny,nu] = size(d);
z = zeros(0,ny);
k = zeros(ny,nu);
zinf = ltipack.getTolerance('infzero',true);
if nu==1
   for i=1:ny
      [zi,gi] = ltipack.sszero(a,b,c(i,:),d(i,:),[],zinf);
      [mz,nz] = size(z);
      nzi = length(zi);
      if i==1,
         z = zi;
      else
         linf = inf;
         z = [[z; linf(ones(max(0,nzi-mz),1),ones(max(nz,1),1))], ...
            [zi;linf(ones(max(0,mz-nzi),1),1)]];
      end
      k(i) = gi;
   end
end