%Program 6.? Animation program for bridge using IVP solver
function tacoma(int,ic,h,p)
%Inputs: int = [a b] time interval,
%ic = [y(1,1) y(1,2) y(1,3) y(1,4)], initialize
%h = stepsize, p = steps per point plotted
%Calls a one-step method such as trapstep.m
%Example usage: tacoma([0 400],[1 0 0.001 0],.04,3)
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;len=6;
 set(gca,'XLim',[-8 8],'YLim',[-8 8], ...
   'XTick',[-8 0 8],'YTick',[-8 0 8], ...
   'Drawmode','fast','Visible','on','NextPlot','add');
cla;              % clear screen
axis square       % make aspect ratio 1 - 1
road =line('color','b','LineStyle','-','LineWidth',5,'erase','xor',...
   'xdata',[],'ydata',[]);
lcable=line('color','r','LineStyle','-','LineWidth',1,'erase','xor',...
   'xdata',[],'ydata',[]);
rcable=line('color','r','LineStyle','-','LineWidth',1,'erase','xor',...
   'xdata',[],'ydata',[]);
for k=1:n
  for i=1:p
    t(i+1) = t(i)+h;
    y(i+1,:) = trapstep(t(i),y(i,:),h);
  end
  y(1,:) = y(p+1,:);t(1)=t(p+1);
z1(k)=y(1,1);z3(k)=y(1,3);tim(k)=t(1);
c=len*cos(y(1,3));s=len*sin(y(1,3));
set(road,'xdata',[c -c],'ydata',[s-y(1,1) -s-y(1,1)])
set(lcable,'xdata',[-c -c],'ydata',[-s-y(1,1) 8])
set(rcable,'xdata',[c c],'ydata',[s-y(1,1) 8])
drawnow; pause(h)
end


function y = trapstep(t,x,h)
%one step of the trapezoid method
z1=ydot(t,x);
g=x+h*z1;
z2=ydot(t+h,g);
y=x+h*(z1+z2)/2;

function ydot=ydot(t,y)
len=6;a=.1;
a1=exp(a*(y(1)-len*sin(y(3))));
a2=exp(a*(y(1)+len*sin(y(3))));
ydot(1) = y(2);
ydot(2) = -0.01*y(2)-0.4*(a1+a2-2)/a+11*sin(3*t);
ydot(3) = y(4);
ydot(4) = -0.01*y(4)+1.2*cos(y(3))*(a1-a2)/(len*a);