%   [I, O, err] = zdsvd(G, Z, epp)
%
%   The ``Zero-Direction-SVD'' function computes the input and
%   output zero directions through the use of Singular Value 
%   Decomposition of G(z) where z is a zero.
%   If singular value number nr (nr=min{ny,nu}, ny-number of outputs, 
%   nu-number of inputs) of G(z) is less than EPP then the 
%   corrosponding input and output singular vectors are stored in
%   columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second smallest singular 
%   value is less than EPP.
%   
%   Inputs:  G   - system matrix.
%            Z   - vector of zeros.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1e-14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to zeros Z(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to zeros Z(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain smaller 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains smaler than EPP.
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   See also: izde, ozde, pdsvd.
%
   function [I, O, err] = zdsvd(G, Z, epp)

   I = []; O = []; err = 0;
   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 3
      epp = 1e-14;
   end

   for i=1:max(size(Z))
      wz  = Z(i)/sqrt(-1);
      Gz = frsp(G, wz);
      [u,s,v] = svd(Gz(1:ny,1:nu));
      nr = min(size(s));
      if s(nr, nr) < epp
         if nr > 1
            if s(nr-1,nr-1) < epp
               disp(['Warning: Second smallest singular value is less than eep: ', num2str(epp),','])
               disp(['         S(',int2str(nr-1),',',int2str(nr-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);       % Rotate elements of vector.
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      else
         disp(['Warning: G(z) is not singular for Z(i): ',num2str(Z(i)), ' for i=',int2str(i),'.'])
         I(:,i) = zeros(nu,1);
	 O(:,i) = zeros(ny,1);
	 err = 1;
      end
   end
end

   
