%Nova representaçao pendulo
clear;clc;close all;

%Definição das Variáveis
% m=0.1; l=0.2; k=0.015/(0.2^2); g=9.81;
g=9.81; m=g^-2; l=g; k=1; 
npont=60;
xp = [-3 3;-5 5];  %definição dos limites dos eixos do gráfico e dos estados
x1min = xp(1,1);  %definição de pontos
x1max = xp(1,2);
x2min = xp(2,1);
x2max = xp(2,2);
deltx1=(x1max-x1min)/npont; % Resolução dos estados
deltx2=(x2max-x2min)/npont;

% for i=1:2
% 
%     x1e(i)=(i-1)*pi;
%     a=[0 1;-g*cos(x1e(i))/l -k/m];  %montando a matriz para calcular os autovalores
%     disp('auto-valor da matriz a para i igual a:')
%     i  %exibir os pontos de equilíbrio
%     eig(a) %calcular os autovalores de A
% end
%resolução das EDO´S

t=[0:0.001:10]; %tempo de resoluçao da EDO
y0=[0.5 1];    %condição inicial para Y
[T,Y]=ode23(@rigid,t,y0);
    
x1=Y(:,1);
x2=Y(:,2);


%plotar V e a solução
figure (1)
q = [-2:deltx1:2];
w = [-5:deltx2:5];
[q,w] = meshgrid(q,w);

%Criar curvas de nivel em 3 dimensões
z=(1-cos(q)).*g/l+w.^2/2;         
contour3 (q,w,z,100,'b');
hold on
plot(Y(1,1),Y(1,2),'r*',Y(:,1),Y(:,2),'r')
xlabel('x1');
hold off
ylabel('x2');
zlabel('Função de Energia');

%Curvas de nível com duas dimensões com trajetória
figure;
contour(q,w,z,50,'b');
hold on
plot(x1,x2,'r')
xlabel('x1');
hold off
ylabel('x2');



%Região de atração com derivada
figure;
v=@(q,w) (1-cos(q))*g/l+w.^2/2;
[q,w] = meshgrid(-0.5:0.001:0.5);
z = v(q,w);
contour(q,w,z,1,'b');
xlim auto
zeroPoints = zeros(1,length(w));
hold on
plot(q(1,:),zeroPoints,'r');
xlabel('x1');
hold off
ylabel('x2');