function temp=zrootinfo(a,ts)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%
% AUTHOR:  Rob McEwen
%
% DATE:    90/5/17
%
% PURPOSE: Display z plane root locations in a comprehensible form.
%
% INPUTS: A:    A matrix 
%         TS:   Sampling period
%
% MODS:
% (1) 93/10/18 rsm.  Log(z) doesn't exist for z <= 0 and z real, which is
%     a legal pole location in the z plane. 
% (2) 94/1/16 rsm.  I found a case where eig(a) gives a much more accurate
%     result than roots(poly(a)).  For this reason, this routine has been
%     modified to accept the matrix itself instead of its characeristic 
%     polynomial.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%
   disp(' ')
   disp('  ---Z Root Information Algorithm---------------------------------')
%
% This lops trailing zeros off the coeff array.  Trailing zeros correspond
% to roots of zero value, for which the log function is undefined.  The
% case for z real and negative has not been handled yet.  93/10/18.
   ncoef = size(a);
   i = ncoef(1,2);
%   while coef(i) = 0, i=i-1; end
%   coef1 = coef(1,1:i);
%   if i < ncoef(1,2),...
%      disp([' There is (are) ' num2str(ncoef(1,2)-i) 'z = 0 roots.']),...
%      disp([' These have been removed.']),...
%   end
%   EVALS   =  ROOT(coef1);
   evalsz   =  eig(a);
   nevalsz  =  size(evalsz);
   nsz      =  nevalsz(1,1);
%
% 94/1/16: remove evals which are zero.
   evals = [0;0];
   j=1;
   for i=1:nsz,...
     if abs(evalsz(i))>1e-7,...
        evals(j)=evalsz(i);...
        j=j+1;...
     end,...
   end
   ns=j-1;

  if ns<nsz,...
  disp([' ' num2str(nsz-ns) 'eigenvalues at z=0 have been removed.']),end

   spl     =  log(evals)/ts;         %s-plane location
   i1=0;
   i2=0;
   i3=0;
%
% split into complex pairs and real roots:
%
   reval = [0;0];
   cmpeval = [0;0];
   for i=1:ns,
     if abs(imag(evals(i)))<=1e-6,
           if abs(real(evals(i)))<=1e-6  i3=i3+1;end,
           i1=i1+1;
           reval(i1)=real(spl(i));
     else
           i2=i2+1;
           cmpeval(i2)=spl(i);
     end
   end
   st1=[' there is (are) ' num2str(i1) ' real evals, ' num2str(i2/2) ...
   ' complex pairs'];
   st2=[' and ' num2str(i3) ' s=0 evals.'];
   disp(st1)
   disp(st2)
%
   if min(abs(evals))>sqrt(2)*1e-7 ...
      disp(' magnitude and phase (deg):'),...
      [abs(evals) asin(  imag(evals)./abs(evals)  )*180/pi],...
   end
   if i2>0  ...
      omn     =  abs(cmpeval);...
      zeta    =  abs(real(cmpeval))./omn;...
      tau     =  -1./real(cmpeval);...
      disp('damping ratio, natural frequency (hz), time constant, 3*t.c.:'),...
      [ zeta omn/2/pi tau 3*tau],...
   end
   if i1>0  ...
      disp('time constant, 3 time constants (seconds):'),...
      tau = [0;0];
      for i=1:i1,...
        if reval(i)>=-1e-5  ...
          tau(i)=0;...
        else,...
          tau(i) = -1/reval(i);...
        end,...
      end,...
      [ tau 3*tau],...
   end
%    tau=-1./.reval;...     fast, but will divde by zero for z=1.
   dummy=0;
   disp('  ---end z root information algorithm------------------------------')
   disp(' ')
return