% A Simulation Comparing Models of the Fall of the World Trade Center Towers
% See http://math.uni-graz.at/keeling/wtc/wtc.html for details

% number of floors:
N = 110;

% crash floor:
C = 100;

% distance between floors in meters:
h = 3.79;

% mass fraction turned to dust:
mu = 0.5;

% floors needed to break next floor:
mcrit = 5;

% gravitational acceleration in meters/seconds^2
g = 9.81;

% height of building:
H = N*h;

% height of crash floor:
X = C*h;

% initially falling mass in units of numbers of floor masses:
minit = N-C+1;

% pcrit is the critical momentum needed to break next floor:
hcrit = 1;
vcrit = sqrt(2*g*hcrit*h);
pcrit = mcrit*vcrit;

% times t(k) at which building has height x(k):
t = zeros(N+1,1);
x = zeros(N+1,1);
t(1) = 0;
x(1) = H;

% velocity and mass of the initially falling object:
mcoll = minit;
vcoll = 0;
% and after falling one floor:
mfall = mcoll;
vfall = -sqrt(vcoll^2 + 2*g*h);

for k=2:C
% time and height of combined object after reaching the next floor:
  t(k) = t(k-1) + (vcoll <= 0.0)*(vcoll + sqrt(vcoll^2 + 2*g*h))/g;
  x(k) = x(k-1) - (vcoll <= 0.0)*h;
% mass and velocity of combined object just after collision:
  mcoll = 1 - mu + mfall;
  vcoll = (pcrit + mfall*vfall)/mcoll;
% mass and velocity of combined object after falling to the next floor:
  mfall = mcoll;
  vfall = -sqrt(vcoll^2 + 2*g*h);
% if the momentum is not sufficient for vcoll<0 then collapse stops:
  if (vcoll > 0)
     t(k+1:end) = t(k);
     x(k+1:end) = x(k);
     break;
  end
end

% the initially falling mass continues in free fall
if (vcoll <= 0)
  for k=(C+1):(N+1)
    t(k) = t(k-1) + (vcoll + sqrt(vcoll^2 + 2*g*h))/g;
    x(k) = x(k-1) - h;
    vcoll = -sqrt(vcoll^2 + 2*g*h);
  end
end

% the free-fall counterpart to the collapse:
y = H - (g/2)*t.^2;
f = sqrt(2*H/g);

disp(strcat('Free fall time=',num2str(f)));
disp(strcat('Pancake Collapse Tim=',num2str(t(end))));

figure(1);
plot(t,N*x/H,'r',t,N*y/H,'b');
axis([t(1) t(N+1) 0 N]);
xlabel('Time (seconds)');
ylabel('Height (floors)');
legend('collapse','free fall');
title('A pancake collapse of a WTC tower is compared with a free fall');
