clc
clear all
close all

%Variaveis do problema
l = 1; %m
k = 1; %N.m/rad = N.m
Pc = k/l; %N

%VetorQ = vetor da coordenada generalizada:
vetorQ = [-90:0.01:90]*pi/180; %rad

%Definicao das cargas nos itens a), b) e c):
%Item a) P_a < Pc
P_a = 0.8*Pc;
%Item b) P_b = Pc
P_b = 1.0*Pc;
%Item b) P_c > Pc
P_c = 1.2*Pc;

%Funcao de Energia Potencial nos itens a) e b) e c):
%Item a) P_a < Pc
PI_a = 1/2*k*vetorQ.^2 - P_a*l*(1 - cos(vetorQ));
%Item b) P_b = Pc
PI_b = 1/2*k*vetorQ.^2 - P_b*l*(1 - cos(vetorQ));
%Item c) P_c > Pc
PI_c = 1/2*k*vetorQ.^2 - P_c*l*(1 - cos(vetorQ));

%Plotando a funcao potencial nos itens a), b) e c):
Fig = figure;
%Item a) P_a < Pc
subplot(3,1,1)
plot(vetorQ.*180/pi, PI_a, 'k', 'linewidth', 1.5)
xlabel('Q[graus]', 'fontsize', 18)
ylabel('\Pi_a[Nm]', 'fontsize', 18)
set(gca,'xlim',[-90,90]) 
%Item b) P_b = Pc
subplot(3,1,2)
plot(vetorQ.*180/pi, PI_b, 'k', 'linewidth', 1.5)
xlabel('Q[graus]', 'fontsize', 18)
ylabel('\Pi_b[Nm]', 'fontsize', 18)
set(gca,'xlim',[-90,90]) 
%Item c) P_c > Pc
subplot(3,1,3)
plot(vetorQ.*180/pi, PI_c, 'k', 'linewidth', 1.5)
xlabel('Q[graus]', 'fontsize', 18)
ylabel('\Pi_c[Nm]', 'fontsize', 18)
set(gca,'xlim',[-90,90]) 

[min,posmin] = min(PI_c)
vetorQ(posmin)*180/pi

##%Output dos resultados
##print('Energias Potenciais', '-dpng', '-r500')