close all;
clear all;

m=1;            % Massa
L=0.2;          % Comprimento
g=9.81;         % Aceleraçao da gravidade
theta0=0;       % Orientaçao inicial
c=0.205;        % Coeficiente de atrito viscoso
k=0.2;          % Coeficiente de flexibilidade

theta0=[0:1e-3:pi];                 % Gama de theta0 (orientaçao de referencia) a testar
res=zeros(length(theta0),2);        % Posiçao dos polos para o sistema sem atrito
res_atr=zeros(length(theta0),2);    % Posiçao dos polos para o sistema com atrito

% Numerador da FT
num=[1/(m*L^2)];

% Varrendo os varios theta0
for i=1:length(theta0)
    % Denominador da FT para o sistema sem atrito
    den=[1 0 g/L*cos(theta0(i))];
    % Denominador da FT para o sistema com atrito
    den_atr=[1 c/(m*L^2) k/(m*L^2)+g/L*cos(theta0(i))];
    
    % Calculo do par de polos para o sistema sem atrito
    res(i,:)=roots(den).';
    % Calculo do par de polos para o sistema com atrito
    res_atr(i,:)=roots(den_atr).';
end

% Figura 1 - Visualizaçao do par de polos para os varios theta0, no sistema
% sem atrito:
figure(1);
% Plot da parte imaginaria em funçao da parte real dos polos
plot(real(res(:,1)),imag(res(:,1)),real(res(:,2)),imag(res(:,2)));
title('Posiçao dos zeros para o pendulo sem atrito - theta0=[0,pi]');
xlabel('Parte Real');
ylabel('Parte Imaginaria');

% Figura 1 - Visualizaçao do par de polos para os varios theta0, no sistema
% com atrito:
figure(2);
% Plot da parte imaginaria em funçao da parte real do segundo polo
plot(real(res_atr(:,1)),imag(res_atr(:,1)),real(res_atr(:,2)),imag(res_atr(:,2)));
title('Posiçao dos zeros para o pendulo com atrito - theta0=[0,pi]');
xlabel('Parte Real');
ylabel('Parte Imaginaria');
