%% prog metodo de euler

clearvars;
close all;
clc;

%% incializacion

%ecuacion diferencial
% y'' + mu*y' + (g/L)*y = 0

%uso de variable auxiliar
% y' = u <------------------------------------> u = y'
% y'' = -mu*y' - (g/L)*y <---------------> u' = -mu*u  - (g/L)*sin(y)

%datos
mu = 0.1;
g = 10.0;
L = 1.0;

%tiempo incial/final
tini =  0;
tfin =  100*pi;

%valores inciales%
y0 = -0.1*pi;     %posición incial
u0 =  0.1*pi;      %velocidad incial u(0) = y'(0)

%incrementos
deltat = 0.0005;

%valores de xn por prelocalización
tn = tini:deltat:tfin;	%valores de tn
Nn = length(tn);        %tamanho de tn

%valores yn por prelocalización
yn = zeros(1,Nn);     %incialmente lleno de ceros

%valores un por prelocalización
un = zeros(1,Nn);     %incialmente lleno de ceros

%condicion incial de desplazamiento
yn(1) = y0;         %condicion incial desplazamiento

%condicion inicial de velocidad
un(1) = u0;         %condicion incial velocidad u(0) = y'(0)


%% solucion numérica
for n = 1:Nn-1
    yn(n+1) = yn(n) + deltat*un(n);
    un(n+1) = un(n) + deltat*( -mu*un(n)  - (g/L)*yn(n) );
end;

%% espacio de fase

Lyini = -1.1*max(abs(yn));
Lyfin =  1.1*max(abs(yn));
dy = 0.1;

Luini = -1.1*max(abs(un));
Lufin =  1.1*max(abs(un));
du = 0.1;

%espacio de estados inciales
[y0,u0] = meshgrid(Lyini:dy:Lyfin,Luini:du:Lufin);

%aceleracion
v0 = -mu*u0  - (g/L)*sin(y0);

%vector
figure(1)
quiver(y0,u0,u0,v0,'r');          %dibuja flechas 2D 
title('espacio de fase');
xlabel('\theta (rad) ');
ylabel('d\theta/dt (rad/s)')
legend('espfase');
daspect([1 1 1]);
axis([ -1.1*max(abs(yn)) 1.1*max(abs(yn)) 1.1*min(un) 1.1*max(un) ]); 
grid on;
box on;

%% resultados

figure(2)
plot(tn,yn,'k')
title('solucion euler desplazamiento');
xlabel('t (s) ')
ylabel('\theta (rad)');
legend(['\theta \Deltat = ' num2str(deltat)])
axis([ tini tfin -1.1*max(abs(yn)) 1.1*max(abs(yn))]); 
grid on;
box on;

figure(3)
plot(tn,un,'k')
title('solucion euler velocidad');
xlabel('t (s)')
ylabel('d\theta/dt (rad/s)')
legend(['d\theta/dt \Deltat = ' num2str(deltat)])
axis([ tini tfin 1.1*min(un) 1.1*max(un) ]); 
grid on;
box on;

figure(4)
plot(yn,un,'k')
title('diagrama de fase');
xlabel('\theta (rad) ');
ylabel('d\theta/dt (rad/s) ');
legend(['\theta vs d\theta/dt \Deltat = ' num2str(deltat)])
axis([ -1.1*max(abs(yn)) 1.1*max(abs(yn)) 1.1*min(un) 1.1*max(un) ]); 
grid on;
box on;

figure(5)
quiver(y0,u0,u0,v0,'r');          %dibuja flechas 2D 
hold on;
plot(yn,un,'k')
title('espacio de fase y diagrama de fase');
xlabel('\theta (rad) $$');
ylabel('d\theta/dt (rad/s) $$');
legend('espfase',['\theta vs d\theta/dt \Deltat = ', num2str(deltat)] );
daspect([1 1 1]);
axis([ -1.1*max(abs(yn)) 1.1*max(abs(yn)) 1.1*min(un) 1.1*max(un) ]); 
grid on;
box on;