function cl_bode_plots
%
% Plots closed-loop transfer functions and responses for several gains
%
%   Outputs: cl_bode_plots.eps       (Frequency magnitude response plots)
%            cl_impulse_response.eps (Impulse responses)
%            cl_step_response.eps    (Step responses)
%
% Actually, this script is meant to be a prototype for a PSTricks plot and
% an example verifying the math.
%

%
% Copyright (c) 2008 by Theodore P. Pavlic
%
% This work is licensed under the Creative Commons
% Attribution-Noncommercial 3.0 United States License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc/3.0/us/
% or send a letter to Creative Commons, 171 Second Street, Suite 300,
% San Francisco, California, 94105, USA.
%

% Upper-case    A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
% Lower-case    a b c d e f g h i j k l m n o p q r s t u v w x y z
% Digits        0 1 2 3 4 5 6 7 8 9
% Exclamation   !           Double quote "          Hash (number) #
% Dollar        $           Percent      %          Ampersand     &
% Acute accent  '           Left paren   (          Right paren   )
% Asterisk      *           Plus         +          Comma         ,
% Minus         -           Point        .          Solidus       /
% Colon         :           Semicolon    ;          Less than     <
% Equals        =           Greater than >          Question mark ?
% At            @           Left bracket [          Backslash     \
% Right bracket ]           Circumflex   ^          Underscore    _
% Grave accent  `           Left brace   {          Vertical bar  |
% Right brace   }           Tilde        ~

%%%%%%%%%%%%%%%%%%%%%%%%%%% Start of Script %%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure(1); clf;
subplot(131); cl_bode_plot(1, -5, -2);
subplot(132); cl_bode_plot(5, -5, -2);
subplot(133); cl_bode_plot(50, -5, -2);

% Uncomment these lines if you want to save the EPS
%
%set( gcf, ...
%    'Position', [0 0 800 600], ...
%    'PaperType', 'usletter', ...
%    'PaperOrientation', 'portrait', ...
%    'PaperPosition', [0.0 0.0 15 4.5] );
%saveas( gcf, 'cl_bode_plots.eps', 'epsc2' );

figure(2); clf;
subplot(131); cl_impulse_response(1, -5, -2);
subplot(132); cl_impulse_response(5, -5, -2);
subplot(133); cl_impulse_response(50, -5, -2);

% Uncomment these lines if you want to save the EPS
%
%set( gcf, ...
%    'Position', [0 0 800 600], ...
%    'PaperType', 'usletter', ...
%    'PaperOrientation', 'portrait', ...
%    'PaperPosition', [0.0 0.0 15 4.5] );
%saveas( gcf, 'cl_impulse_response.eps', 'epsc2' );

figure(3); clf;
subplot(131); cl_step_response(1, -5, -2); 
subplot(132); cl_step_response(5, -5, -2); 
subplot(133); cl_step_response(50, -5, -2); 

% Uncomment these lines if you want to save the EPS
%
%set( gcf, ...
%    'Position', [0 0 800 600], ...
%    'PaperType', 'usletter', ...
%    'PaperOrientation', 'portrait', ...
%    'PaperPosition', [0.0 0.0 15 4.5] );
%saveas( gcf, 'cl_step_response.eps', 'epsc2' );


    function [sigma, dampfreq, zeta] = second_order_closed_loop( K, a, b )

        % Open-loop parameters are K (gain), a and b (poles)

        % Closed-loop poles
        p = roots( [1 -(a+b) (1+K)*a*b ] );

        % Damping and damped frequency
        sigma = -real(p(1));
        dampfreq = imag(p(1));

        % Damping ratio (sigma/(natural frequency))
        % zeta = sigma(1)/abs(p(1));
        zeta = sigma(1)./sqrt( sigma(1).^2 + dampfreq(1).^2 );

    end

    function cl_bode_plot(K, a, b)
        
        % Frequency range, on a log scale (i.e., range of exponents of 10)
        logw = linspace(0,2,1000);

        % Open-loop parameters are K (gain), a and b (poles)

        % Find real and imaginary parts of CL poles and damping ratio
        [sigma, dampfreq, zeta] = second_order_closed_loop( K, a, b );
        
        % Plot the closed-loop bode
        % (note: The elegant MATLAB way to do this would be to use complex
        %        numbers, the abs function, and a loglog plot. However, 
        %        I'm trying to generate a nice prototype for pst-plot. 
        %        Plus, doing things "raw" like this might be a good
        %        example.)
        plot( logw, ...
              20*log10( K*a*b./sqrt( ((1+K)*a*b - (10.^logw).^2).^2 + (a+b)^2*(10.^logw).^2 ) ), ...
              '-', ...
              logw, 0*logw, '--' );
        ylim( [-20 15] );
        grid on;
        
        if( length(zeta) > 0 )
            legstr = [ '$-\sigma = ', num2str(-sigma(1)), '$', ...
                       '\\$j\omega_d = j', num2str(dampfreq(1)), '$', ...
                       '\\$\zeta = ', num2str(zeta(1), 3) , '$' ];
        else
            
            if( a < b )
                legstr = [ '$-\sigma_a = ', num2str( -sigma(1) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(2) ), '$' ];
            else
                legstr = [ '$-\sigma_a = ', num2str( -sigma(2) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(1) ), '$' ];
            end
                       
        end
        
        legstr = [ legstr, ...
                   '\\$e_{SS} = ', num2str(1/(1+K)), '$' ];

        title( [ '$K = ' num2str(K) '$, ' ...
                 '$a = ' num2str(a) '$, ' ...
                 '$b = ' num2str(b) '$' ], ...
               'Interpreter', 'latex', 'FontSize', 14 );

        xlabel( '$\log_{10}(\omega/(1$\,rad$/$s$))$', ...
                'Interpreter', 'latex', 'FontSize', 12 );
        ylabel( '$|G_{CL}|$ (dB)', ...
                'Interpreter', 'latex', 'FontSize', 12 );

        lfig = legend( '', 'Location', 'NorthWest' );
        set( lfig, 'Interpreter', 'latex', ...
                   'FontSize', 12, ...
                   'String', [ '\shortstack{' legstr '}' ] );
            
    end

    function cl_impulse_response(K, a, b, Tmax)
        
        if( nargin < 4 )
            Tmax = 2;
        end
        
        % Time range
        t = linspace(0, Tmax, 2048);
        
        % Open-loop parameters are K (gain), a and b (poles)

        % Find real and imaginary parts of CL poles and damping ratio
        [sigma, dampfreq, zeta] = second_order_closed_loop( K, a, b );
        
        % Plot the closed-loop bode
        % (note: The elegant MATLAB way to do this would be to use the
        %        toolboxes (e.g., the impulse function).)
        plot( t, ...
              K/(1+K) * ...
              (dampfreq^2+sigma^2) / dampfreq * ...
              exp(-sigma*t) .* ...
              sin(dampfreq*t), ...
              '-', ...
              t, 0*t, '--' );
        ylim( [-10 20] );
        grid on;
        
        if( length(zeta) > 0 )
            legstr = [ '$-\sigma = ', num2str(-sigma(1)), '$', ...
                       '\\$j\omega_d = j', num2str(dampfreq(1)), '$', ...
                       '\\$\zeta = ', num2str(zeta(1), 3) , '$' ];
        else
            
            if( a < b )
                legstr = [ '$-\sigma_a = ', num2str( -sigma(1) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(2) ), '$' ];
            else
                legstr = [ '$-\sigma_a = ', num2str( -sigma(2) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(1) ), '$' ];
            end
                       
        end
        
        legstr = [ legstr, ...
                   '\\$e_{SS} = ', num2str(1/(1+K)), '$' ];

        title( [ '$K = ' num2str(K) '$, ' ...
                 '$a = ' num2str(a) '$, ' ...
                 '$b = ' num2str(b) '$' ], ...
               'Interpreter', 'latex', 'FontSize', 14 );

        xlabel( 'Time (s)', ...
                'Interpreter', 'latex', 'FontSize', 12 );
        ylabel( 'Impulse response (e.g., velocity profile)', ...
                'Interpreter', 'latex', 'FontSize', 12 );

        lfig = legend( '', 'Location', 'NorthEast' );
        set( lfig, 'Interpreter', 'latex', ...
                   'FontSize', 12, ...
                   'String', [ '\shortstack{' legstr '}' ] );
            
    end

    function cl_step_response(K, a, b, Tmax)
        
        if( nargin < 4 )
            Tmax = 2;
        end
        
        % Time range
        t = linspace(0, Tmax, 2048);
        
        % Open-loop parameters are K (gain), a and b (poles)

        % Find real and imaginary parts of CL poles and damping ratio
        [sigma, dampfreq, zeta] = second_order_closed_loop( K, a, b );
        
        % Plot the closed-loop bode
        % (note: The elegant MATLAB way to do this would be to use the
        %        toolboxes (e.g., the impulse function).)
        plot( t, ...
              K/(1+K) * ...
                 ( 1 - ...
                       ( cos(dampfreq*t) - ...
                         (sigma/dampfreq)*sin(dampfreq*t) ) .* ...
                       exp(-sigma*t) ), '-', ...
              t, K/(1+K)*ones(size(t)), '--', ...
              t, ones(size(t)), '-.' );
        ylim( [0 2] );
        grid on;
        
        if( length(zeta) > 0 )
            legstr = [ '$-\sigma = ', num2str(-sigma(1)), '$', ...
                       '\\$j\omega_d = j', num2str(dampfreq(1)), '$', ...
                       '\\$\zeta = ', num2str(zeta(1), 3) , '$' ];
        else
            
            if( a < b )
                legstr = [ '$-\sigma_a = ', num2str( -sigma(1) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(2) ), '$' ];
            else
                legstr = [ '$-\sigma_a = ', num2str( -sigma(2) ), '$', ...
                           '\\$-\sigma_b = ', num2str( -sigma(1) ), '$' ];
            end
                       
        end
        
        legstr = [ legstr, ...
                   '\\$e_{SS} = ', num2str(1/(1+K)), '$' ];

        title( [ '$K = ' num2str(K) '$, ' ...
                 '$a = ' num2str(a) '$, ' ...
                 '$b = ' num2str(b) '$' ], ...
               'Interpreter', 'latex', 'FontSize', 14 );

        xlabel( 'Time (s)', ...
                'Interpreter', 'latex', 'FontSize', 12 );
        ylabel( 'Impulse response (e.g., velocity profile)', ...
                'Interpreter', 'latex', 'FontSize', 12 );

        lfig = legend( '', 'Location', 'NorthEast' );
        set( lfig, 'Interpreter', 'latex', ...
                   'FontSize', 12, ...
                   'String', [ '\shortstack{' legstr '}' ] );
            
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%% End of Script %%%%%%%%%%%%%%%%%%%%%%%%%%%