%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% forcepend.m
% Animation program for forced pendulum 
% using IVP solver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function forcepend(int, ic, h, p, d, force)

% Inputs: int = [a b] time interval,

% ic = [y(1,1) y(1,2)], initialize

% h = stepsize, p = steps per point plotted

% Calls a one-step method such as trapstep.m

% Example usage: pend([0 10],[pi/2 0],.05,1)

clf                       % clear figure window

a=int(1);b=int(2);n=ceil((b-a)/(h*p)); % plot n points in total

y(1,:)=ic;                % enter initial conds in y

t(1)=a;xmin=10;

set(gca,'XLim',[-1.2 1.2],'YLim',[-1.2 1.2], ...

  'XTick',[-1 0 1],'YTick',[-1 0 1], ...

  'Drawmode','fast','Visible','on','NextPlot','add');

cla;            % clear screen

plot(0,0,'ks')  % pivot where rod attached

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',[]);

for k=1:n

  for i=1:p

    t(i+1) = t(i)+h;

    y(i+1,:) = forcetrapstep(t(i),y(i,:),h, d, force);

  end

  y(1,:) = y(p+1,:);t(1)=t(p+1);

  xbob = cos(y(1,1)-pi/2); ybob = sin(y(1,1)-pi/2);xmin=min(xbob,xmin);

  xrod = [0 xbob]; yrod = [0 ybob];

  set(rod,'xdata',xrod,'ydata',yrod)

  set(bob,'xdata',xbob,'ydata',ybob)

  drawnow; pause(h)

end

xmin