clear all;
clc;

%% Create the parameters

bqr = 1; % Assign the value of q*b^2/r

a = -20:1e-3:20;
acl = zeros(size(a));
P = zeros(size(a));
acl_den = zeros(size(a));
i = 1;

for aa = a
    
    Lambda(i) = bqr/(1-aa^2);
    if abs(aa)>1
        acl(i) = aa/(1-((1-aa^2)/2)*(-Lambda(i)+1+sqrt((1-Lambda(i))^2+4*Lambda(i)/(1-aa^2))));
    elseif abs(aa)<1
        acl(i) = aa/(1+((1-aa^2)/2)*(-Lambda(i)+1+sqrt((1-Lambda(i))^2+4*Lambda(i)/(1-aa^2))));
    end

    i = i+1;
    
end

% Plot the closed-loop system gain
figure(1);
plot(a,acl, 'LineWidth', 2);
set(gca, 'FontSize',15,'fontname','times','xlim',[min(a), max(a)],'ylim',[min(acl)-0.3,max(acl)+0.3]);
grid on;
xlabel('$$\mathbf{a}$$','Interpreter','Latex');
ylabel('$$\mathbf{a^{cl}}$$','Interpreter','Latex');
title('$$\textbf{Closed-Loop System Gain}$$','Interpreter','Latex')