clear; close all
%
%  Two masses connected by spring & damper
%  Dropped to ground
%
m1 = 100;
m2 = 10;
k = 4000;
b = 25;
duration = 15;
dt = duration/800;
t = 0:dt:duration;
height = 0.5;
g = 9.81;
%
%  Initial conditions
%
x0 = height*[1;1];
xdot0 = [0;0];
y0 = [x0; xdot0];
%
%  Don't allow x2 < 0
%
options = odeset('NonNegative',2);
%[t,y] = ode45(@twomassfun,t,y0,options,[m1,m2,k,b]);  %% Revised Mar. 3
[t,y] = ode45(@twomassfun,t,y0,options,m1,m2,k,b);  %% Revised Mar. 3
%
%  Plot mass motions vs. time
%
plot(t,1000*y(:,1));
hold on
plot(t,1000*y(:,2),'r');
legend('Upper mass','Lower mass');
xlabel('Time, sec');
ylabel('Position, mm');
%
%  Draw system, initial condition
%
figure
size1 = 0.1;
size2 = 0.05;
separation = 1;
hb1 = 0;
c1 = [0 height+separation+size1/2];
hb1 = drawbox(c1,size1,size1,'b');
ydata1 = get(hb1,'YData');
hold on
c2 = [0 height+size2/2];
hb2 = drawbox(c2,size1,size2,'b');
ydata2 = get(hb2,'YData');
hk = plot(-[1 1]*size2/3,[c1(2)-size1/2 c2(2)+size2/2],'b','LineWidth',4);
hb = plot(+[1 1]*size2/3,[c1(2)-size1/2 c2(2)+size2/2],'g','LineWidth',4);
ydatak = get(hk,'YData');
set(gca,'YLim',[-height 2*height+separation]);
set(gca,'XLim',height*[-1 1]);
set(gca,'XTick',[],'YTick',[]);
plot([-1 1]*height/2,[0 0],'k','LineWidth',4);
%
%  Animate motion
%
x1 = y(:,1);
x2 = y(:,2);
nothing = input('Hit ENTER: ');
nt = length(t);
ht = xlabel(sprintf('%8.3f sec',0));
for i=1:nt
    set(hb1,'YData',y(i,1)-height+ydata1);
    set(hb2,'YData',y(i,2)-height+ydata2)
    pause(0.01);
    set(hk,'YData',[separation+x1(i) size2+x2(i)]);
    set(hb,'YData',[separation+x1(i) size2+x2(i)]);
    set(ht,'String',sprintf('%8.3f sec',t(i)));
end;