function animDoublePend(t,y)
%Function to animate the pendulum.  Input: t and y
%  from the numerical solution to the pendulum equations.
%  Length is hard-coded to be 1.
[m,n]=size(y);
n2=length(t);

%Error check on sizes (y should be "tall")
if m<n
    y=y';
end
[m,n]=size(y);
if m~=n2
    error('Time values and Position do not match');
end

plot(0,0,'ks');           % pivot where rod attached
axis square               % make aspect ratio 1 - 1
bob1 = line('color','k','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod1 = line('color','b','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);
bob2 = line('color','r','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod2 = line('color','g','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);
for k=1:m-1
    theta1=y(k+1,1)-pi/2; theta2=y(k+1,3)-pi/2;
    xbob1 = cos(theta1);       ybob1 = sin(theta1);
    xbob2 = xbob1+cos(theta2); ybob2=ybob1+sin(theta2);
    xrod1 = [0 xbob1]; yrod1 = [0 ybob1];
    xrod2 = [xbob1 xbob2]; yrod2 = [ybob1 ybob2];
    set(rod1,'xdata',xrod1,'ydata',yrod1)
    set(bob1,'xdata',xbob1,'ydata',ybob1)
    set(rod2,'xdata',xrod2,'ydata',yrod2);
    set(bob2,'xdata',xbob2,'ydata',ybob2);
    axis([-2.2 2.2 -2.2 2.2]);
    drawnow; pause(0.01)
end