function [prop,eigA]=systprop(A,B,C,D)
% SYSTPROP - computes the following properties of a linear system:
%
%          1. time constant                             tau
%          2. natural frequency of undamped system      w0
%          3. eigenfrequency of the system              wn
%          4. period                                    P
%          5. damping factor                            zeta
%          6. percentage overshoot                      PO
%          7. peak-time                                 Tpeak
%          8. settling time to a 2% band                Tset
%          9. halve-time                                Thalve
%
% Conclusions are displayed on screen and stored in a user-specified
% diary-file (by default SYSTPROP.LOG). This includes information about
% the stability of the system and possible oscillatory behaviour, and
% a short summary of the meaning of some system properties. 
%
% The default location in which the diary-file will be stored is the 
% FDC data directory; if SYSTPROP is not able to locate that directory, 
% it will default to the current working directory. You can change the 
% destination directory and filename in the save dialog. Creation of 
% the diary-file can be skipped altogether by selecting 'Cancel' or 
% by clicking the Close button in the save dialog. 
%
% [X,Y] = SYSTPROP(A,B,C,D) or [X,Y] = SYSTPROP(A,B,C) returns the 
%    system properties in the matrix X and the eigenvalues of A in the 
%    columnvector Y. In addition, this call returns information about 
%    the controllability and observability of the system.
%
%    The colums of the matrix X correspond with the system properties 
%    1 through 9 in the table above: 
%
%      X = [tau, w0, wn, P, zeta, PO, Tpeak, Tset, Thalve]
%
%    If the elements of the D matrix are all zero, it isn't necessary to
%    enter the D matrix; in that case you may use [X,Y] = SYSTPROP(A,B,C).
%
% [X,Y] = SYSTPROP(A) returns the system properties in the matrix X and 
%    the eigenvalues of A in the columvector Y; see the definitions above.
%
% [X,Y] = SYSTPROP(num,den) evaluates the equivalent transfer function
%    representation, again returning the matrix X and vector Y. It also
%    checks whether there are zeros in the right half of the s-plane, 
%    to identify minimum-phase and non-minimum-phase systems.
%
% X = SYSTPROP(Y) returns the system properties in the matrix X if Y 
%    contains the eigenvalues of a system in a columnvector. If Y is a 
%    scalar, it will be treated as a single eigenvalue, instead of a
%    1x1 A-matrix. To analyze a system consisting of a 1x1 A-matrix, 
%    use the call [X,Y] = SYSTPROP(A,0,0) instead.
%
% Note: this function requires the Control System Toolbox! In addition, 
% be aware that SYSTPROP uses the routine NUM2STR2 for output formatting.
%
% See also DAMP, CTRB, OBSV.


% ---------------------------------------------------------------------
% REFERENCE: J. van de Vegte, 'Feedback Control Systems', Prentice Hall
%            International Editions, London, 2nd edition, 1990.
% ---------------------------------------------------------------------


% Disable Matlab warning messages and clear the command-window, to prevent the 
% screen output from being distorted. This function was originally written for 
% Matlab 3.5, and although it has been updated since then, it may cause inadver-
% tent warning messages due to slightly obsolete coding. 
% ------------------------------------------------------------------------------
warning off
clc

nargs=nargin;

% If number of inputs is zero: return help information only. Else, compute
% the system properties and store the results in a diary-file.
% ------------------------------------------------------------------------
if nargs == 0
    helpwin systprop
else

    % Diary-file specification. The while-loop gives the user the possibility to
    % return to filename-specification when the cancel-button of the uiputfile GUI
    % is pressed inadvertantly. The variable usediary determines whether or not
    % a diary-file will be used. By default it is 'on'. This variable can be 
    % changed only by editing systprop.m. 
    % ----------------------------------------------------------------------------
    usediary = 'on';
    ok = 'No';
    while not(strcmp(ok,'Yes'))

        % Let the user specify a suitable name and location for the SYSTPROP diary-file.
        % The default file-name is SYSPROP.LOG; the default location is the FDC data-
        % directory, or, if that directory can't be found, the current working directory.
        % -------------------------------------------------------------------------------
        defaultdir = feval('datadir');
        workingdir = pwd;

        if exist(defaultdir,'dir')
            [filename,pathname] = uiputfile([defaultdir filesep 'systprop.log'], ...
                                            'Select name and location of diary-file:');
        else
            [filename,pathname] = uiputfile([workingdir filesep 'systprop.log'], ...
                                            'Select name and location of diary-file:');
        end

        % If cancel-button was clicked, ask whether or not to omit creating the diary-file.
        % If the user doesn't want a diary-file, change the variable usediary to 'off'.
        % ---------------------------------------------------------------------------------
        if not(ischar(filename) | ischar(pathname))
            ok= questdlg(['No diary-file has been specified. Continue without diary-file?'], ...
                          'WARNING', ...
                          'Yes','No','No');
            if strcmp(ok,'Yes')
                usediary = 'off';
            end
        else
            ok = 'Yes';
        end

    end

    % If required, activate specified diary-file. If file already exists, delete
    % it (a question whether or not to replace an existing file will be created
    % by uiputfile itself, so we don't need to take additional precautions here).
    % ---------------------------------------------------------------------------
    if strcmp(usediary,'on')
        location = [pathname filename];
        if exist(location,'file')
            delete(location);
        end
        diary(location);
    end


    % Now determine which actions to take, depending on the number of input arguments
    % -------------------------------------------------------------------------------
    if nargs >= 3              % Three or four inputs: treat inputs as system matrices
                               % of a linear state-space system [A,B,C,D]

        % If the D matrix is not specified, create D-matrix as a matrix of zeros
        % ----------------------------------------------------------------------
        if nargs == 3;
            [m1,n1] = size(B);
            [m2,n2] = size(C);
            D = zeros(m2,n1);
        end

        % Display system information as specified in SYSTPROP inputs
        % ----------------------------------------------------------
        disp('Linear state-space system with system matrices A, B, C, D:');
        A,B,C,D
        disp(' ');

        % Display error message if the dimensions of the system matrices don't match
        % --------------------------------------------------------------------------
        error(abcdchk(A,B,C,D));

        % Display information about controllability and observability of the system
        % -------------------------------------------------------------------------
        Co = ctrb(A,B);
        Ob = obsv(A,C);
        rankA = rank(A);

        disp(' ');
        disp(' ');
        if rank(Co) == rankA & rank(Ob) == rankA
            disp('System is controllable and observable.');
        elseif rank(Co) == rankA
            disp('System is controllable.');
        elseif rank(Ob) == rankA
            disp('System is observable.');
        else
            disp('System is neither controllable nor observable');
        end;

        % Determine eigenvalues of the system matrix A
        % --------------------------------------------
        eigA = eig(A);

    elseif nargs == 1          % Only one input: treat input as system matrix A or as a vector
                               % with eigenvalues Y. One scalar input will always be treated as
                               % a single eigenvalue instead of a 1x1 matrix A. It is still
                               % possible to input a 1x1 matrix A by defining the matrices B,
                               % C, and D too, e.g.: SYSTPROP(A,0,0,0) is a valid call.
        [m,n] = size(A);

        if n==1                % Treat input A as a vector with eigenvalues.

            % Display system information as specified in SYSTPROP input
            % ---------------------------------------------------------
            Y = A;
            disp('Vector with eigenvalues Y:');
            Y

            % Store vector with eigenvalues in eigA
            % -------------------------------------
            eigA = A;

        elseif m == n          % Treat input A as state matrix.

            % Display system information as specified in SYSTPROP input
            % ---------------------------------------------------------
            disp('System matrix A:');
            A

            % Determine eigenvalues of the system matrix A
            % --------------------------------------------
            eigA = eig(A);

        else
            error('Error:  Incorrect input-type for SYSTPROP');
        end

    elseif nargs == 2          % Two inputs: treat inputs A and B as numerator and
                               % denominator of a transfer function. In this case, A
                               % is a vector with coefficients of the numerator and B
                               % is a vector with coefficients of the denominator.

        % Display system information as specified in SYSTPROP inputs
        % ----------------------------------------------------------
        num = A;
        den = B;
        tf(num,den)

        % Determine poles of the system
        % -----------------------------
        eigA = roots(den);
        eigB = roots(num);
        if max(sign(eigB)) == 1
            disp('Zeros found in right half-plane => non-minimum phase system.')
        elseif max(sign(eigB)) == -1
            disp('All zeros in left half-plane => minimum phase system.')
        end
        disp(' '); 

    else                       % More than four inputs: not correct.
        error('Error: too many inputs for SYSTPROP')
    end


    % Find additional system properties by analyzing the eigenvalues of the 
    % system. First, determine the number of system modes.
    % ---------------------------------------------------------------------
    [m,n] = size(eigA);
    
    if m == 1
       disp('There is one system mode:');
    else
       disp(['There are ' num2str(m) ' system modes:']);
    end
    disp(' ');

    % Next, determine the system properties for each individual mode
    % --------------------------------------------------------------
    for i = 1:1:m
        ra = real(eigA(i,1));                             % Real parts of the eigenvalues
        ia = imag(eigA(i,1));                             % Imaginary parts of the eigenvalues

        prop(i,9) =  log(0.5)/ra;                         % halve time [s]
        prop(i,1) = -1/ra;                                % time constant [s]
        prop(i,2) =  abs(eigA(i,1));                      % undamped natural frequency [rad/s]
        prop(i,5) = -ra/prop(i,2);                        % damping coefficient [-]
        prop(i,8) =  4*prop(i,1);                         % settling time [s]

        % If the eigenvalues have a non-zero imaginary part, determine the properties of the
        % oscillatory signals. For non-oscillatory behavior, the eigenfrequency equals zero,
        % the period is infinite, as is the peak time (which equals pi/eigenfrequency = pi/0).
        % ------------------------------------------------------------------------------------
        if ia ~= 0
            prop(i,3)= prop(i,2)*sqrt(1-prop(i,5)^2);     % eigenfrequency [rad]
            prop(i,4)= abs(2*pi/prop(i,3));               % period [s]
            prop(i,7)= pi/prop(i,3);                      % peak time [s]
        else
            prop(i,3)= 0;                                 % eigenfrequency [rad]
            prop(i,4)= inf;                               % 'period' [s]
            prop(i,7)= inf;                               % 'peak time' [s]
        end

        % Determine the overshoot percentage
        % ----------------------------------
        if prop(i,5) > 0 & prop(i,5) < 1;                 % damped oscillatory
            prop(i,6) = 100*exp(-pi*prop(i,5)/sqrt(1-prop(i,5)));
        elseif prop(i,5) > 1;                            % non-oscillatory and stable
            prop(i,6) = 0;
        else
            prop(i,6) = inf;                              % unstable
        end 

        % Analyze the results (categorize the characteristic responses for the different
        % eigenvalues into corresponding 'system modes' depending on their eigenvalues).
        %-------------------------------------------------------------------------------
        if ia == 0 & ra < 0
            mode ='stable, non-oscillatory';
        elseif ia == 0 & ra > 0
            mode ='unstable, non-oscillatory';
        elseif ia ~= 0 & ra < 0
            mode ='stable, oscillatory';
        elseif (ia ~= 0 & ra > 0) | (prop(i,5) < 0 & prop(i,5) > -1)
            mode ='unstable, oscillatory';
        elseif prop(i,5) == 1 | ra == 0
            mode ='critically damped, indifferent';
        elseif prop(i,5) < -1
            mode ='unstable, non-oscillatory';
        elseif prop(i,5) > 1
            mode ='non-oscillatory';
        else
            mode ='unknown (!!)';
        end;

        % Display system modes conclusions
        % --------------------------------
        if ia == 0
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra) ]);
        elseif ia < 0
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra), ...
                  ' ' num2str(ia) 'i']);
        else
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra), ...
                  ' +' num2str(ia) 'i']);
        end

    end

    NewPage(usediary); % local function, see below

    % Display results. Note: the routine NUM2STR2 is called instead of the standard Matlab 
    % routine NUM2STR in order to obtain better screen formatting. For 1, 2, or 3 system
    % modes, the page formatting is tuned for optimal readibility. If there are more than
    % 3 system modes, a short-hand notation will be used.
    %-------------------------------------------------------------------------------------
    if m == 1

        disp('                          Mode 1')
        disp('                          ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8) '   [s]']);
        disp(['Undamped nat. freq. :   ' num2str2(prop(1,2),8) '   [rad/s]']);
        disp(['Eigenfrequency      :   ' num2str2(prop(1,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8) '   [s]']);

    elseif m == 2

        disp('                          Mode 1       Mode 2')
        disp('                          ======       ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8), ...
              '     ' num2str2(prop(2,1),8) '   [s]']);
        disp(['Natural frequency   :   ' num2str2(prop(1,2),8), ...
              '     ' num2str2(prop(2,2),8) '   [rad/s]']);
        disp(['Eigenfrequency      :   ' num2str2(prop(1,3),8), ...
              '     ' num2str2(prop(2,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8), ...
              '     ' num2str2(prop(2,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8), ...
              '     ' num2str2(prop(2,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8), ...
              '     ' num2str2(prop(2,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8), ...
              '     ' num2str2(prop(2,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8), ...
              '     ' num2str2(prop(2,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8), ...
              '     ' num2str2(prop(2,9),8) '   [s]']);

    elseif m == 3

        disp('                          Mode 1       Mode 2       Mode 3')
        disp('                          ======       ======       ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8), ...
              '     ' num2str2(prop(2,1),8) '     ' num2str2(prop(3,1),8) '   [s]']);
        disp(['Natural frequency   :   ' num2str2(prop(1,2),8), ...
              '     ' num2str2(prop(2,2),8) '     ' num2str2(prop(3,2),8) '   [rad/s]']);
        disp(['Eigen frequency     :   ' num2str2(prop(1,3),8), ...
              '     ' num2str2(prop(2,3),8) '     ' num2str2(prop(3,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8), ...
              '     ' num2str2(prop(2,4),8) '     ' num2str2(prop(3,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8), ...
              '     ' num2str2(prop(2,5),8) '     ' num2str2(prop(3,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8), ...
              '     ' num2str2(prop(2,6),8) '     ' num2str2(prop(3,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8), ...
              '     ' num2str2(prop(2,7),8) '     ' num2str2(prop(3,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8), ...
              '     ' num2str2(prop(2,8),8) '     ' num2str2(prop(3,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8), ...
              '     ' num2str2(prop(2,9),8) '     ' num2str2(prop(3,9),8) '   [s]']);

    elseif m > 3

        disp('           tau [s]  w0 [rad/s]  wn [rad/s]       P [s]    zeta [-]');
        for i = 1:1:m
            disp(['Mode ' int2str(i) ':   ' num2str2(prop(i,1),8), ...
                  '    ' num2str2(prop(i,2),8) '    ' num2str2(prop(i,3),8), ...
                  '    ' num2str2(prop(i,4),8) '    ' num2str2(prop(i,5),8)]);
        end
        disp(' ');
        disp('            PO [%]   Tpeak [s]    Tset [s]  Thalve [s]');
        for i=1:1:m
            disp(['Mode ' int2str(i) ':   ' num2str2(prop(i,6),8), ...
                  '    ' num2str2(prop(i,7),8) '    ' num2str2(prop(i,8),8), ...
                  '    ' num2str2(prop(i,9),8)]);
        end
        disp(' ');
        disp('tau   : time constant                      PO    : percentage overshoot');
        disp('w0    : natural freq. of undamped system   Tpeak : peak-time');
        disp('wn    : eigenfrequency of the system       Tset  : settling-time (2%)');
        disp('P     : period                             Thalve: halve-time');
        disp('zeta  : damping factor');
        
    else               % Not 1, 2, 3, or more system modes? Then something must be wrong. 
        error('Invalid number of system modes!');
    end

    disp(' ');

    NewPage(usediary); % local function, see below

    disp('Notes:')
    disp('======')
    disp('  1. The system is stable (i.e. the transient solution will decay to zero)')
    disp('     if all poles are in the left half plane of the s-plane')
    disp('  2. To avoid excessive overshoot and unduly oscillator behaviour, the')
    disp('     damping ratio must be adequate, thus the angle phi (the angle between the')
    disp('     real-axis of the s-plane and the line through the origin and the pole)')
    disp('     must not be too close to zero.')
    disp('  3. The time constant and settling time can be reduced (response speed')
    disp('     increased) by increasing the negative real part of the poles.')
    disp('  4. Undamped frequency = distance from pole to origin; moving the poles out')
    disp('     radially (i.e. with constant damping factor) increases the speed of');
    disp('     response, thus reducing settling time, peak time and rise time, while the')
    disp('     percentage overshoot remains constant.')
    disp('  5. Natural frequency = resonant frequency or damped natural frequency. This')
    disp('     is the frequency of transient oscillations, which equals the imaginary')
    disp('     part of the pole positions.')
    disp('  6. Peak time and rise time are reduced by increasing the imaginary part')
    disp('     of the pole locations.')
    disp(' ');
    disp('Refer to source-code of SYSTPROP.M for more details about the computations.')
    disp('See also: J. van de Vegte, ''Feedback Control Systems'', Prentice Hall')
    disp('International Editions, London, 2nd edition, 1990.')
    
    NewPage(usediary); % local function, see below
    
    disp('Ready.')
    disp(' ')

    if strcmp(usediary,'on')
        diary off
        disp('The results have been stored in the diary-file:');
        disp(location);
    end
end

% Enable Matlab warning messages again
% ------------------------------------
warning on


% ====================
% BEGIN LOCAL FUNCTION
% ====================

    function NewPage(usediary)
    % Small local function, turns off diary function when required (that is: when
    % usediary equals 'on'), waits for user to press a key, clears the screen and
    % turns diary function back on when required. This prevents the 'press a key to
    % continue' messages to be included in the SYSTPROP logfile. 

    disp(' ');
    if strcmp(usediary,'on')
        diary off
    end
    disp('<<< Press a key to continue >>>');
    pause
    clc
    if strcmp(usediary,'on')
        diary on
    end

% ==================
% END LOCAL FUNCTION
% ==================


%-----------------------------------------------------------------------
% This program is based on a tool written by E.A. van der Zwan in 1993.
%
% The Flight Dynamics and Control Toolbox version 1.4.0. 
% (c) Copyright Marc Rauw, 1994-2005. Licensed under the Open Software 
% License version 2.1; see COPYING.TXT and LICENSE.TXT for more details.
% Last revision of this program: July 18, 2005.
%-----------------------------------------------------------------------
