% Exemplo 16.17 - H(s) = (s^2+10s+100)/((s+1)(s+100)), xi = 0.5 e wn=10
clear all
close all
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotando o gráfico do módulo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w = logspace(-2, 4);
xi = 0.5;
wn=10;
for i = 1:length(w)
    MdB(i) = mag2db(sqrt((wn^2-w(i)^2)^2+(2*xi*wn*w(i))^2))+mag2db(1/sqrt(w(i)^2+1))+mag2db(1/sqrt(w(i)^2+100^2)); 
end
figure
semilogx(w,MdB,'r','LineWidth',2);
axis([10^(-2) 10^4 -40 40])
grid on
hold on

% Plotando assíntotas 
% wc=1:
w = [10^(-2) 10^0]; 
M = [0 0];
plot(w,M,'b:','LineWidth',2);

w = [10^(0) 10^4]; 
M = [-20*log10(w(1)) -20*log10(w(2))];
plot(w,M,'b:','LineWidth',2);

% wc=100:
w = [10^(-2) 10^2]; 
M = [0 0];
plot(w,M,'b:','LineWidth',2);

w = [10^(2) 10^4]; 
M = [40-20*log10(w(1)) 40-20*log10(w(2))];
plot(w,M,'b:','LineWidth',2);

% wc=10:
w = [10^(-2) 10^1]; 
M = [0 0];
plot(w,M,'b:','LineWidth',2);

w = [10^(1) 10^4]; 
M = [-40+40*log10(w(1)) -40+40*log10(w(2))];
plot(w,M,'b:','LineWidth',2);

% H = H1H2H3
w = [10^(-2) 10^0];
M = [0 0];
plot(w,M,'b','LineWidth',2);

w = [10^(0) 10^1];
M = [-20*log10(w(1)) -20*log10(w(2))];
plot(w,M,'b','LineWidth',2);

w = [10^(1) 10^2];
M = [-20*log10(w(1))-40+40*log10(w(1)) -20*log10(w(2))-40+40*log10(w(2))];
plot(w,M,'b','LineWidth',2);

w = [10^(2) 10^4];
M = [-20*log10(w(1))-40+40*log10(w(1))+40-20*log10(w(1)) -20*log10(w(2))-40+40*log10(w(2))+40-20*log10(w(2))];
plot(w,M,'b','LineWidth',2);


title('Diagrama assintótico de módulo')
xlabel('\omega')
ylabel('M_{dB}(\omega)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotando o gráfico da fase
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w = logspace(-2, 4);
for i = 1:length(w)
    if w(i)<=wn
        Phi(i) = (-atan(w(i))-atan(w(i)/100)+atan(2*xi*wn*w(i)/(wn^2-w(i)^2)))*180/pi;
    elseif w(i)>wn
        Phi(i) = (-atan(w(i))-atan(w(i)/100)+pi+atan(-2*wn*xi*w(i)/(w(i)^2-wn^2)))*180/pi;
    end
end
figure
semilogx(w,Phi,'r','LineWidth',2);
axis([10^(-2) 10^4 -200 200])
grid on
hold on


% Plotando assíntotas de fase
% (wc=1) 
w = [10^(-2) 10^(-1)]; 
PA = [0 0];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(-1) 10^(1)]; 
PA = [0 -90];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(1) 10^(4)]; 
PA = [-90 -90];         
plot(w,PA,'b:','LineWidth',2);

% (wc=100) 
w = [10^(-2) 10^(1)]; 
PA = [0 0];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(1) 10^(3)]; 
PA = [0 -90];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(3) 10^(4)]; 
PA = [-90 -90];         
plot(w,PA,'b:','LineWidth',2);

% (wc=10) (raízes complexas - mudança abrupta em w=wc):
w = [10^(-2) 10^(1)]; 
PA = [0 0];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(1) 10^(1)]; 
PA = [0 180];         
plot(w,PA,'b:','LineWidth',2);

w = [10^(1) 10^(4)]; 
PA = [180 180];         
plot(w,PA,'b:','LineWidth',2);

% H=H1H2H3
w = [10^(-2) 10^(-1)]; 
PA = [0 0];         
plot(w,PA,'b','LineWidth',2);

w = [10^(-1) 10^(1)]; 
PA = [0 -90];         
plot(w,PA,'b','LineWidth',2);

w = [10^(1) 10^(1)]; 
PA = [-90 90];         
plot(w,PA,'b','LineWidth',2);

w = [10^(1) 10^(3)]; 
PA = [90 0];         
plot(w,PA,'b','LineWidth',2);

w = [10^(3) 10^(4)]; 
PA = [0 0];         
plot(w,PA,'b','LineWidth',2);

title('Diagrama assintótico de fase')
xlabel('\omega')
ylabel('\Phi(\omega)')
