% Exemplo 16.14 - Polos e zeros reais
clear all
close all
clc

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

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

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

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

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

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

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

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

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

w = [10^(2) 10^3]; 
M = [-40 -40];         
plot(w,M,'b','LineWidth',2);

w = [10^3 10^5];   
M = [20-20*log10(w(1)) 20-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, 5);
for i = 1:length(w)
    Phi(i) = (-atan(w(i))+atan(w(i)/100)-atan(w(i)/1000))*180/pi; 
end
figure
semilogx(w,Phi,'k','LineWidth',2);
axis([10^(-2) 10^5 -100 100])
grid on
hold on

% Plotando assíntotas considerando H=H1*H2*H3
% H1 (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^(5)]; 
PA = [-90 -90];         
plot(w,PA,'b:','LineWidth',2);

% H2 (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^(5)]; 
PA = [90 90];         
plot(w,PA,'b:','LineWidth',2);

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

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

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

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

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

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

w = [10^2 10^3];   
P = [-45 -45];     
plot(w,P,'b','LineWidth',2);

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

w = [10^4 10^5];   
P = [-90 -90];     
plot(w,P,'b','LineWidth',2);


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