% Program 6.3 Animation program for pendulum
% Inputs: time interval inter,
%  initial values ic = [y(1,1) y(1,2)], number of steps n
% Calls a one-step method such as trapstep.m
% Example usage: pend([0 10],[pi/2 0],200)
function exctish(inter,ic,n)
h=(inter(2)-inter(1))/n;  % plot n points in total
set(gca,'xlim',[-2.7 2.7],'ylim',[-2.7 2.7], ...
  'XTick',[-2.5 0 2.5],'YTick',[-2.5 0 2.5], ...
  'Drawmode','fast','Visible','on','NextPlot','add');
cla;         
axis square               % make aspect ratio 1 - 1
bob = line('color','r','Marker','.','markersize',40,...
    'erase','xor','xdata',[],'ydata',[]);
rod = line('color','b','LineStyle','-','LineWidth',3,...
    'erase','xor','xdata',[],'ydata',[]);
[t,y] = ode45(@ydot,inter, ic);
v=VideoWriter('exctsh_stbl','MPEG-4');open(v)
for k=1:n
  frame=getframe;writeVideo(v,frame);
  xbob = 2.5*sin(y(k+1,1)); ybob = -2.5*cos(y(k+1,1));
  %xbob
  xrod = [0 xbob]; yrod = [0 ybob];
  set(rod,'xdata',xrod,'ydata',yrod)
  set(bob,'xdata',xbob,'ydata',ybob)
  xlabel(['A = ' num2str(18.356)])
  drawnow; pause(h)
end
close(v)

function z=ydot(t,y)
g=9.81;length=2.5;
z = zeros(2,1);
z(1) = y(2);
z(2) = -((g/length)+ 18.356*cos(2*pi*t))*sin(y(1)) - 0.1*y(2) ;