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.4;          % Coeficiente de flexibilidade

theta0=[0:pi/10: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: sistema sem atrito/flexibilidade nas juntas
    den=[1 0 g/L*cos(theta0(i))];
    % Denominador da FT: sistema com atrito e flexibilidade nas juntas
    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:
%subplot(2,1,1)
figure(1);
plot(real(res(:,1)),imag(res(:,1)),'bx', real(res(:,2)),imag(res(:,2)),'bx',...
     [-10 10],[0 0],'k', [0 0],[-10 10],'k') 
title('Sistema sem atrito e flexibilidade nas juntas');
xlabel('Sigma');
ylabel('jw');

% Figura 1 - Visualizaçao do par de polos para os varios theta0, no sistema
% com atrito e flexibilidade nas juntas:
%subplot(2,1,2)
figure(2)
plot(real(res_atr(:,1)),imag(res_atr(:,1)),'bx', real(res_atr(:,2)),imag(res_atr(:,2)),'bx',...
    [-10 10],[0 0],'k', [0 0],[-10 10],'k') 
title('Sistema com atrito e flexibilidade nas juntas');
xlabel('Sigma');
ylabel('jw');

