clear all
close all
clc

% Para H(s)=1/(s^2+s+1)=w0^2/(s^2+2*xi*w0*s+w0^2)
% tem-se que:
% wo^2=1
w0 = 1;
% 2*xi*w0=1
xi = 1/(2*w0);

w = 0:0.001:5; % vetor de frequências para plotar o gráfico

for i=1:length(w)
    % Calcula o módulo de H(jw) pela expressão obtida:
    Mw(i) = w0^2/sqrt((w0^2-w(i)^2)^2 + 4*xi^2*w0^2*w(i)^2);
    % Calcula a fase de H(jw) pela expressão obtida:
    Phiw(i) = -atan(w(i))*180/pi;      
end

% Desenha os gráficos:
figure(1)
plot(w,Mw)
xlabel('\omega')
ylabel('|H(j\omega)|')
grid on

Mwr = max(Mw); % Calcula o ponto de máximo do módulo pelo gráfico
ind = find(Mw==Mwr); % encontra a posição no vetor Mw desse ponto de máximo
wr = w(ind);% encontra o valor de w correspondente ao ponto de máximo

% Imprime na tela o módulo máximo e a frequência equivalente
fprintf('Analisando o gráfico, tem-se:\n')
fprintf('M(wr) = %1.4f\n',Mwr)
fprintf('wr = %1.4f\n',wr)
% Imprime na tela o resultado de acordo com a fórmula encontrada na letra b
fprintf('\n')
fprintf('b)\n')
fprintf('M(wr) = %1.4f\n',1/(2*xi*sqrt(1-xi^2)))
fprintf('wr = %1.4f\n',w0*sqrt(1-2*xi^2))

hold on % mantém a figura aberta
plot(wr,Mwr,'ro') % plota sobre a figura o ponto de máximo com um círculo vermelho


    