## Copyright (C) 2017 Mark Bronsfeld
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see .
## -*- texinfo -*-
## @deftypefn {Function File} damp(@var{sys})
## @deftypefnx {Function File} {[@var{Wn}, @var{zeta}] =} damp(@var{sys})
## @deftypefnx {Function File} {[@var{Wn}, @var{zeta}, @var{P}] =} damp(@var{sys})
## Calculate natural frequencies, damping ratios and poles.
##
## If no output is specified, display overview table containing poles,
## magnitude (if @var{sys} is a discrete-time model), damping ratios, natural
## frequencies and time constants.
##
## @strong{Inputs}
## @table @var
## @item sys
## @acronym{LTI} model.
## @end table
##
## @strong{Outputs}
## @table @var
## @item Wn
## Natural frequencies of each pole of @var{sys} (in increasing order).
## The frequency unit is rad/s (radians per second).
## If @var{sys} is a discrete-time model with specified sample time, @var{Wn}
## contains the natural frequencies of the equivalent continuous-time poles
## (see Algorithms). If @var{sys} has an unspecified sample time
## (@var{tsam =} -1), @var{tsam =} 1 is used to calculate @var{Wn}.
## @item zeta
## Damping ratios of each pole of @var{sys} (in the same order as @var{Wn}).
## If @var{sys} is a discrete-time model with specified sample time,
## @var{zeta} contains the damping ratios of the equivalent continuous-time
## poles (see Algorithms). If @var{sys} has an unspecified sample time
## (@var{tsam =} -1), @var{tsam =} 1 is used to calculate @var{zeta}.
## @item P
## Poles of @var{sys} (in the same order as @var{Wn}).
## @end table
##
## @strong{Algorithm}@*
## @table @var
## @item Pole
## Poles s (or z for discrete-time models) are calculated via @command{pole}
## and resorted in order of increasing natural frequency.
## @item Equivalent continuous-time pole
## s = log (z) / sys.tsam (discrete-time models only)
## @item Magnitude
## mag = abs (z) (discrete-time models only)
## @item Natural Frequency
## Wn = abs (s)
## @item Damping ratio
## zeta = -cos (arg (s))
## @item Time constant
## tau = 1 / (Wn * zeta)
## @end table
##
## @seealso{dsort, eig, esort, pole, pzmap, zero}
## @end deftypefn
## Author: Mark Bronsfeld
## Created: January 2017
## Version: 0.1
function [Wn_out, zeta, P] = damp (sys)
if (nargin == 1) # damp (sys)
if (! (isa (sys, "lti") || issquare (sys)))
error ("damp: argument must be an LTI system");
endif
else
print_usage ();
endif
P = pole (sys); # Poles
## Distinguish between system/state matrices, continuous- and
## discrete-time models
if ((! (isa (sys, "lti")) && issquare (sys)) || (isct (sys)))
s = P;
elseif (isdt (sys))
if (sys.tsam == -1) # If sample time is unspecified...
s = log (P); # ...assume 1 second: log (P) / 1
else
s = log (P) ./ sys.tsam;
endif
mag = abs (P); # Magnitude
endif
Wn = abs (s); # Frequencies (rad / seconds)
## Sort all vectors in order of increasing natural frequency
[Wn, ndx] = sort (Wn);
P = P (ndx);
s = s (ndx);
zeta = -cos (arg (s)); # Damping
tau = 1 ./ (Wn .* zeta); # Time constant (seconds)
## Suppress "ans" output if no output specified (only assign "Wn_out"
## if any output specified)
if (nargout > 0)
Wn_out = Wn;
## Display overview table when no output specified
elseif (nargout == 0)
## Type conversion and formatting to exponential format
P = num2str (P, '%1.2e');
zeta = num2str (zeta, '%1.2e');
Wn = num2str (Wn, '%1.2e');
tau = num2str (tau, '%1.2e');
## Construct columns of overview table
Pole = [['Pole'; ' '; ' ']; P];
Damping = [['Damping'; ' '; ' ']; zeta];
Frequency = [['Frequency'; '(rad/seconds)'; ' ']; Wn];
TimeConstant = [['Time Constant'; '(seconds)'; ' ']; tau];
## Construct overview table - distinguish between
## system/state matrices, continuous- and discrete-time models
if ((! (isa (sys, "lti")) && issquare (sys)) || (isct (sys)))
## Overview table
overview_table = [repmat(' ',rows(Pole),3) Pole repmat(' ',rows(Pole),4) Damping repmat(' ',rows(Pole),4) Frequency repmat(' ',rows(Pole),4) TimeConstant];
elseif (isdt (sys))
mag = mag (ndx); # Sort vector in order of increasing natural frequency
mag = num2str (mag, '%1.2e'); # Type conversion and formatting to exponential format
Magnitude = [['Magnitude'; ' '; ' ']; mag]; # Construct additional column of overview table
## Overview table
overview_table = [repmat(' ',rows(Pole),3) Pole repmat(' ',rows(Pole),4) Magnitude repmat(' ',rows(Pole),4) Damping repmat(' ',rows(Pole),4) Frequency repmat(' ',rows(Pole),4) TimeConstant];
endif
disp (overview_table); # Display overview table
endif
endfunction
## Test system/state matrix
%!test
%! A = [-1, 0, 0;
%! 0, -2, 0;
%! 0, 0, -3 ];
%!
%! Wn_exp = [1;
%! 2;
%! 3 ];
%!
%! zeta_exp = [1;
%! 1;
%! 1 ];
%!
%! P_exp = [-1;
%! -2;
%! -3 ];
%!
%! [Wn_obs, zeta_obs, P_obs] = damp (A);
%!
%! assert (Wn_obs, Wn_exp, 0);
%! assert (zeta_obs, zeta_exp, 0);
%! assert (P_obs, P_exp, 0);
## Test continuous-time model
%!test
%! H = tf ([2, 5, 1], [1, 2, 3]);
%!
%! Wn_exp = [1.7321;
%! 1.7321 ];
%!
%! zeta_exp = [0.5774;
%! 0.5774 ];
%!
%! P_exp = [-1.0000 + 1.4142i;
%! -1.0000 - 1.4142i ];
%!
%! [Wn_obs, zeta_obs, P_obs] = damp (H);
%!
%! assert (Wn_obs, Wn_exp, 1e-4);
%! assert (zeta_obs, zeta_exp, 1e-4);
%! assert (P_obs, P_exp, 1e-4);
## Test discrete-time model
%!test
%! H = tf ([5, 3, 1], [1, 6, 4, 4], 0.01);
%!
%! Wn_exp = [193.4924;
%! 193.4924;
%! 356.5264 ];
%!
%! zeta_exp = [ 0.0774;
%! 0.0774;
%! -0.4728 ];
%!
%! P_exp = [-0.3020 + 0.8063i;
%! -0.3020 - 0.8063i;
%! -5.3961 + 0.0000i ];
%!
%! [Wn_obs, zeta_obs, P_obs] = damp (H);
%!
%! assert (Wn_obs, Wn_exp, 1e-4);
%! assert (zeta_obs, zeta_exp, 1e-4);
%! assert (P_obs, P_exp, 1e-4);