%% Analysis of Feedback Control Systems
%
%
% Table of Contents
% Introduction
%     Gain and Phase Margins
%     Issues with the tf( ) Command
%         Excerpt from the online documentation regarding tf( )
%         Proposed Solution
% The Complementary Sensitivity Function
%     Static Gain and Pole-Zero Analysis
%         Static Gain
%         Poles & Zeros
%     Frequency Response Analysis
%     More Examples
% The Sensitivity Function
%% Introduction
% Consider the following open-loop transfer functions

clear
addpath(genpath('BodeDiagram/'))

% let's define the transfer functions
s=tf('s');

L1s = 5/(1+s)^3;
L2s = 0.5*((1+5*s)^2/s)*(1/(1+s)^3);
L3s = 10*(s/(1+s/100))*(1/(1+s)^3);

% and assign a label (a "name") to each one
set(L1s,"Name", "L_1");
set(L2s,"Name", "L_2");
set(L3s,"Name", "L_3");

%
%% Gain and Phase Margins
% How to compute the gain and the phase margins for each of the given transfer functions?

help margin

% Thus, we have for L1s

[Gm1, Pm1, OMpi1, OMc1] = margin(L1s);
fprintf(['*** Transfer function L1s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L1s ***\n\r'], ...
         Gm1, OMpi1, Pm1, OMc1);

% Then for L2s:

[Gm2, Pm2, OMpi2, OMc2] = margin(L2s);
fprintf(['*** Transfer function L2s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L2s ***\n\r'], ...
         Gm2, OMpi2, Pm2, OMc2);

%
% And finally for L3s

[Gm3, Pm3, OMpi3, OMc3] = margin(L3s);
fprintf(['*** Transfer function L3s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n', ...
         'Ang. freq. OMpi [rad/s]: %.4f \n\r', ...
         'Phase margin [deg]: %.4f \r\n',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L3s ***\n\r'], ...
         Gm3, OMpi3, Pm3, OMc3);

%
%
%% Issues with the tf( ) Command
%
% Pay attention! If you use the following (correct) syntax, you will encounter possible numerical problems. In fact, the tf() command is NOT able to handle simplifications correctly.

F1sALT = L1s/(1+L1s)
zpk(F1sALT)

%
%% Excerpt from the online documentation regarding tf( )
%
% Limitations
% Transfer function models are ill-suited for numerical computations. Once created, convert them to state-space form before combining them with other models or performing model transformations. You can then convert the resulting models back to transfer function form for inspection purposes.
%
%% Proposed Solution
% To avoid these problems, it is preferable to use the feedback() command to obtain the closed-loop transfer function from the open-loop one.

F1s = feedback(L1s, 1, -1)

%
%% The Complementary Sensitivity Function
%
% Given the transfer functions L_1(s), \, L_2(s) and L_3(s), the corresponding closed-loop functions F(s) are
% F(s) = \frac{L(s)}{1+L(s)}

F1s = feedback(L1s, 1, -1)
set(F1s,"Name", "F_1");

F2s = feedback(L2s, 1, -1)
set(F2s,"Name", "F_2");

F3s = feedback(L3s, 1, -1)
set(F3s,"Name", "F_3");

%
%% Static Gain and Pole-Zero Analysis
% Remember:

zero(L1s)
pole(L1s)
zero(L2s)
pole(L2s)
pole(L3s)
zero(L3s)

%% Static Gain


dcF10 = dcgain(F1s)

dcF20 = dcgain(F2s)

dcF30 = dcgain(F3s)

%
%% Poles & Zeros

zpk(L1s)

zpk(F1s)


figure('Units','centimeters', 'Position',[0.1, 0.1, 24, 20]);
hPZP = pzplot(F1s, 'r', L1s, 'g'); 
grid on
% what are the properties of a pzplot object?
get(hPZP)
% let's set the legend visible and select the "best possible" position
hPZP.LegendVisible = "on";
hPZP.LegendLocation = "best";

% let's configure color and size for the zeros & poles marker fo the first
% transfer function
get(hPZP.Responses(1))

% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP.Responses(1).MarkerSize = 12;
hPZP.Responses(1).LineWidth = 2;
hPZP.Responses(1).Color = 'r';

% set font and size of axes labels
get(hPZP.AxesStyle(1))

% see 'doc pzplot' & 'doc pzoptions'
hPZP.AxesStyle(1).FontSize = 14;
hPZP.AxesStyle(1).GridColor = [0,0,0];

hPZP.XLabel(1)

hPZP.XLabel(1).FontSize = 14;
hPZP.YLabel(1).FontSize = 14;

% let's configure color and marker size for  zeros & poles of the second
% transfer function
hPZP.Responses(2).MarkerSize = 12;
hPZP.Responses(2).LineWidth = 2;
hPZP.Responses(2).Color = 'g';


%

zpk(L2s)

zpk(F2s)


figure('Units','centimeters', 'Position',[0.1, 0.1, 24, 20]);
hPZP2 = pzplot(F2s, 'r', L2s, 'g'); 
grid on

% let's set the legend visible and select the "best possible" position
hPZP2.LegendVisible = "on";
hPZP2.LegendLocation = "best";

% let's configure color and size for the zeros & poles marker fo the first
% transfer function

% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP2.Responses(1).MarkerSize = 12;
hPZP2.Responses(1).LineWidth = 2;
hPZP2.Responses(1).Color = 'r';
% let's configure color and marker size for  zeros & poles of the second
% transfer function
hPZP2.Responses(2).MarkerSize = 14;
hPZP2.Responses(2).LineWidth = 2;
hPZP2.Responses(2).Color = 'g';
% see 'doc pzplot' & 'doc pzoptions'
hPZP2.AxesStyle(1).FontSize = 14;
hPZP2.AxesStyle(1).GridColor = [0,0,0];

hPZP2.XLabel(1).FontSize = 14;
hPZP2.YLabel(1).FontSize = 14;

%

zpk(L3s)

zpk(F3s)


figure('Units','centimeters', 'Position',[0.1, 0.1, 24, 20]);
hPZP3 = pzplot(F3s, 'r', L3s, 'g'); 
grid on

% let's set the legend visible and select the "best possible" position
hPZP3.LegendVisible = "on";
hPZP3.LegendLocation = "best";

% let's configure color and size for the zeros & poles marker fo the first
% transfer function

% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP3.Responses(1).MarkerSize = 12;
hPZP3.Responses(1).LineWidth = 2.5;
hPZP3.Responses(1).Color = 'r';
% let's configure color and marker size for  zeros & poles of the second
% transfer function
hPZP3.Responses(2).MarkerSize = 16;
hPZP3.Responses(2).LineWidth = 1.5;
hPZP3.Responses(2).Color = 'g';
% see 'doc pzplot' & 'doc pzoptions'
hPZP3.AxesStyle(1).FontSize = 14;
hPZP3.AxesStyle(1).GridColor = [0,0,0];

hPZP3.XLabel(1).FontSize = 14;
hPZP3.YLabel(1).FontSize = 14;

%
%% Frequency Response Analysis
%
% Consider
%
% L_4(s) = \mu_L \, \frac{\left(1+\frac{s}{50}\right)}{ \left(1+\frac{s}{2}\right)\,\left(1+0.8\, s +\frac{s^2}{4}\right)}}
% with \mu_L = 1.2,\, 5.45
% Let's compare the frequency response of the open-loop transfer function L_4(s) and that of the closed-loop transfer function F_4(s).
%

muL4default = 1.20;
muL4 = 5.45; % try muL4 = 5.45
L4s = muL4*(1+s/50)/((1+s/2)*(1+0.8*s+s^2/4));

F4s = feedback(L4s,1,-1);

figure;margin(L4s)

figure;
hPZP4 = pzplot(F4s);
grid on

% let's set the legend visible and select the "best possible" position
hPZP4.LegendVisible = "on";
hPZP4.LegendLocation = "best";

% let's configure color and size for the zeros & poles marker fo the first
% transfer function

% PZPlot properties:
% https://it.mathworks.com/help/ident/ref/controllib.chart.pzplot-properties.html
% https://it.mathworks.com/matlabcentral/answers/2173352-change-marker-colour-of-pzplot
%
hPZP4.Responses(1).MarkerSize = 12;
hPZP4.Responses(1).LineWidth = 2.5;
hPZP4.Responses(1).Color = 'r';
% see 'doc pzplot' & 'doc pzoptions'
hPZP4.AxesStyle(1).FontSize = 14;
hPZP4.AxesStyle(1).GridColor = [0,0,0];

hPZP4.XLabel(1).FontSize = 14;
hPZP4.YLabel(1).FontSize = 14;


pole(F4s)

hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L4s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(F4s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{4}(s)$', 'actual $L_{4}(s)$', ...
             'asympt. $F_{4}(s)$', 'actual $F_{4}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{4}(s)$', 'actual $L_{4}(s)$', ...
             'asympt. $F_{4}(s)$', 'actual $F_{4}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%%
%
%% More Examples
%

[Gm1, Pm1, OMpi1, OMc1] = margin(L1s);
fprintf(['*** Transfer function L1s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L1s ***\n\r'], ...
         Gm1, OMpi1, Pm1, OMc1);

hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L1s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(F1s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $F_{1}(s)$', 'actual $F_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $F_{1}(s)$', 'actual $F_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%

[Gm2, Pm2, OMpi2, OMc2] = margin(L2s);
fprintf(['*** Transfer function L2s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L2s ***\n\r'], ...
         Gm2, OMpi2, Pm2, OMc2);

hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L2s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(F2s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{2}(s)$', 'actual $L_{2}(s)$', ...
             'asympt. $F_{2}(s)$', 'actual $F_{2}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{2}(s)$', 'actual $L_{2}(s)$', ...
             'asympt. $F_{2}(s)$', 'actual $F_{2}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%

[Gm3, Pm3, OMpi3, OMc3] = margin(L3s);
fprintf(['*** Transfer function L3s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n', ...
         'Ang. freq. OMpi [rad/s]: %.4f \n\r', ...
         'Phase margin [deg]: %.4f \r\n',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L3s ***\n\r'], ...
         Gm3, OMpi3, Pm3, OMc3);


hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L3s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(F3s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{3}(s)$', 'actual $L_{3}(s)$', ...
             'asympt. $F_{3}(s)$', 'actual $F_{3}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{3}(s)$', 'actual $L_{3}(s)$', ...
             'asympt. $F_{3}(s)$', 'actual $F_{3}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%

hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L1s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(F1s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $F_{1}(s)$', 'actual $F_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $F_{1}(s)$', 'actual $F_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%
%
%
%% The Sensitivity Function
%

[Gm1, Pm1, OMpi1, OMc1] = margin(L1s);
fprintf(['*** Transfer function L1s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L1s ***\n\r'], ...
         Gm1, OMpi1, Pm1, OMc1);

S1s = tf(1)-F1s
zpk(S1s)
pole(S1s)
zero(S1s)


hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L1s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(S1s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $S_{1}(s)$', 'actual $S_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{1}(s)$', 'actual $L_{1}(s)$', ...
             'asympt. $S_{1}(s)$', 'actual $S_{1}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%
%

[Gm2, Pm2, OMpi2, OMc2] = margin(L2s);
fprintf(['*** Transfer function L1s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L1s ***\n\r'], ...
         Gm2, OMpi2, Pm2, OMc2);

S2s = tf(1)-F2s
zpk(S2s)
pole(S2s)
zero(S2s)


hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L2s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(S2s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{2}(s)$', 'actual $L_{2}(s)$', ...
             'asympt. $S_{2}(s)$', 'actual $S_{2}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2, 'asympt. $L_{2}(s)$', 'actual $L_{2}(s)$', ...
             'asympt. $S_{2}(s)$', 'actual $S_{2}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%

[Gm3, Pm3, OMpi3, OMc3] = margin(L3s);
fprintf(['*** Transfer function L1s ***\n\r\n\r' ...
         'Gain Margin [dB]: %.4f \n\r', ...
         'Ang. freq. OMpi [rad/s]: %.4f \r\n', ...
         'Phase margin [deg]: %.4f \n\r',...
         'Crit. Ang. freq. OMc [rad/s]: %.4f \n\r\n\r',...
         '*** Transfer function L1s ***\n\r'], ...
         Gm3, OMpi3, Pm3, OMc3);

S3s = tf(1)-F3s
zpk(S3s)
pole(S3s)
zero(S3s)


hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);

omVALS = logspace(-1,2, 1e4); 
% We want to evaluate the frequency response at 10000 
% angular pulsation values between 10^-1 and 10^2 rad/s.
% 
% ** Note **: 
%    The greater the number of angular pulsation values 
%    used to evaluate the frequency response, the better
%    the accuracy of the graphs of both the asymptotic 
%    and actual diagrams.

% let's create the diagrams: we need to store the handler of the figure 
% refer to the help of drawBodediagrams( )

asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
           0.9290 0.6940 0.1250; ...
           0.4940 0.1840 0.5560; ...
           0.4660 0.6740 0.1880; ...
           0.3010 0.7450 0.9330; ...
           0.8500 0.3250 0.0980; ...
           0.6350 0.0780 0.1840]; % some different colors 
                                  % for the actual diagrams



drawBodediagrams(L3s, omVALS, asBcolor, 2.5, '-', ...
               Bcolors(1, :), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(S3s, omVALS, Bcolors(2, :), 2.5, '-', ...
               Bcolors(4, :), 1.5, '-', hf);

legend(hax1, 'asympt. $L_{3}(s)$', 'actual $L_{3}(s)$', ...
             'asympt. $S_{3}(s)$', 'actual $S_{3}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');
legend(hax2,  'asympt. $L_{3}(s)$', 'actual $L_{3}(s)$', ...
             'asympt. $S_{3}(s)$', 'actual $S_{3}(s)$', ...
             'Location', 'best', 'Interpreter', 'latex');

%
%
%
%
%
%