function [wnout,z,r] = damp(a,Ts)
%DAMP  Natural frequency and damping of linear system dynamics.
%
%    [Wn,Z] = DAMP(SYS) returns vectors Wn and Z containing the natural
%    frequencies and damping factors of the linear system SYS. For
%    discrete-time models, the equivalent s-plane natural frequency and
%    damping ratio of an eigenvalue lambda are:
%
%       Wn = abs(log(lambda))/Ts ,   Z = -cos(angle(log(lambda))) .
%
%    If the sample time, Ts, is undefined, the software assumes Ts = 1.
%
%    [Wn,Z,P] = DAMP(SYS) also returns the poles P of SYS.
%
%    When invoked without left-hand arguments, DAMP prints the poles with
%    their natural frequencies, damping factors, and time constants in a
%    tabular format on the screen. The poles are sorted by increasing
%    frequency. Both Wn and P are expressed in the reciprocal of the time
%    units of SYS. For example, 1/minute if SYS.TimeUnit = 'minutes'.
%
%    See also POLE, ESORT, DSORT, PZMAP, ZERO, LTI, DYNAMICSYSTEM.

% Old help
%warning(['This calling syntax for ' mfilename ' will not be supported in the future.'])
%DAMP   Natural frequency and damping factor for continuous systems.
%   [Wn,Z] = DAMP(A) returns vectors Wn and Z containing the
%   natural frequencies and damping factors of A.   The variable A
%   can be in one of several formats:
%
%       1) If A is square, it is assumed to be the state-space
%          "A" matrix.
%       2) If A is a row vector, it is assumed to be a vector of
%          the polynomial coefficients from a transfer function.
%       3) If A is a column vector, it is assumed to contain
%          root locations.
%
%   When invoked without left hand arguments DAMP prints the 
%   eigenvalues with their natural frequency and damping factor in a
%   tabular format on the screen.
%
%   See also EIG and DDAMP.

%   J.N. Little 10-11-85
%   Revised 3-12-87 JNL
%   Revised 7-23-90 Clay M. Thompson
%   Copyright 1986-2011 The MathWorks, Inc.
if nargin>0 && ~isnumeric(a)
   error(message('Control:general:NotSupportedModelsofClass','damp',class(a)))
end

sa = size(a);
if sa(2)>1
   if sa(1)==sa(2)
      % first input is an A matrix
      r = eig(a);
   elseif sa(1)==1
      % first input is a NUM vector
      r = roots(a);
   else
      error(message('Control:analysis:damp2'))
   end
else
   r = a;
end

% Sample time argument (undocumented)
if nargin==1
   Ts = 0;
elseif ~isscalar(Ts)
   error(message('Control:analysis:damp3'))
end
Ts = abs(Ts);

% Compute Wn and Z
z = NaN(size(r));
if Ts>0
   % Discrete: compute equivalent S-plane poles
   s = log(r)/Ts;
   z(r==0) = +1;     % r=0 stable
   z(r==1 | isinf(r)) = -1;  % unstable
else
   s = r;
   z(r==0 | isinf(r)) = -1;  % unstable by convention
end
wn = abs(s);
idx = find(wn>0 & wn<Inf);
z(idx) = -real(s(idx))./wn(idx);

% Output
if nargout==0      
	% Print results on the screen. 
	ltipack.printdamp(r,wn,z,Ts)
else
	wnout = wn; 
end


