%% run a simulation of a double pendulum
t = 0:.01:15;
x0 = [pi/2, 0, 0, 0];
[t,x] = ode45('DoublePendulum', t, x0);

%% plot states over time
figure
plot(t,x)

%% run an animation
figure('position',[62 855 965 333])
subplot(1,2,1)
plot(t,x(:,1:2))
legend({'a1','a2'})
hold on
s1ylim = get(gca,'ylim');
tline = plot([0 0],s1ylim,'r-','linewidth',1);
subplot(1,2,2)
[h0,h1,h2] = a2h(x(1,1),x(1,2));
l1 = plot([h0(1) h1(1)],[h0(2) h1(2)],'b-');
hold on
l2 = plot([h1(1) h2(1)],[h1(2) h2(2)],'b-');
axis([-2.2 2.2 -2.2 2.2]);
axis equal
p0 = plot(h0(1),h0(2),'k.','markersize',15);
p1 = plot(h1(1),h1(2),'b.','markersize',15);
p2 = plot(h2(1),h2(2),'r.','markersize',15);
tp = title(sprintf('%5.1f',t(1)));
skip = 3;
pause
for i=1:skip:length(t)
subplot(1,2,1)
set(tline, 'Xdata', [t(i) t(i)]);
set(tline, 'Ydata', s1ylim); subplot(1,2,2)
[h0,h1,h2] = a2h(x(i,1), x(i,2));
set(l1, 'XData', [h0(1) h1(1)]);
set(l1, 'YData', [h0(2) h1(2)]);
set(l2, 'XData', [h1(1) h2(1)]);
set(l2, 'YData', [h1(2) h2(2)]);
set(p1, 'XData', h1(1));
set(p1, 'YData', h1(2));
set(p2, 'XData', h2(1));
set(p2, 'YData', h2(2));
set(tp, 'String', sprintf('%5.1f', t(i)));
axis([-2.2 2.2 -2.2 2.2]);
drawnow;
end