% pendulo simple sin rozamiento
set(0,'DefaultLineLineWidth',3);
clf;
% tiempo
t=[0:1/(24*4):20];
% angulo inicial: probar pequenos (pi/6) y grandes (pi/1.1)
a0=pi/3;
% longitud
l0=10;
% velocidad angular inicial
v0=0;
% masa
m=2;
% coef. rozamiento: 0: no rozamiento (periodico)
r=0;
% gravedad
g=9.81;

% 2 Ley de Newton angular para el pendulo simple
% a es el angulo, w la velocidad angular
% a' = w
% w' = -g/l*sin(a) - r*l/m*w
f=@(t, X)...
    [X(2);...                % esto es a'=w
     -g/l0*sin(X(1)) - r/m*l0*X(2)];  % esto es v'=-g/l*sin(a)-r/m*l*w

% solucion
Y = heunvector(f,t,[a0;v0]);
Yeuler = eulervector(f,t,[a0;v0]);
% angulo en azul
% velocidad angular en rojo
plot(t, Y);


pause;
l=[0:.01:l0];
T=[0:.01:2*pi];
clf
xlim([-l0 l0])
ylim([0.8 2*l0+1])
%hold on
for k=1:floor(length(Y(1,:))/4);
    alfa=Y(1,4*k);
    alfaE=Yeuler(1,4*k);
    %hold off;
    pause(1/24)
    plot([0 l0*cos(3*pi/2-a0)],[l0+1 l0+1+l0*sin(3*pi/2-a0)],'g--','LineWidth',1);
    hold on;
    plot([0 l0*cos(3*pi/2-alfa)],[l0+1 (l0)+1+l0*sin(3*pi/2-alfa)],'b');
    plot([0 l0*cos(3*pi/2-alfaE)],[l0+1 (l0)+1+l0*sin(3*pi/2-alfaE)],'k','LineWidth',1);
    rectangle('Position',[l0*cos(3*pi/2-alfa)-l0/20,l0+1+l0*sin(3*pi/2-alfa)-l0/20 l0/10 l0/10],...
    'Curvature',[1 1],'FaceColor',[0.9 0 0]);
    hold off;
    axis equal
    xlim([-l0 l0])
    ylim([0.8 2*l0+1])    %pause;
    %        hold on;
    %
end
    
    
    