% Parameters for Earthquake Problem/Newmark integration

dt = 0.02;
nsteps = 600 ;
gamma = 0.50 ;
beta = 0.25 ;
alpha1 = 1.2;
alpha2 = 0.0015;
time = i*dt;
% Form Mass, Stiffness and Load Matrices

mass = 1500*[1,0,0,0;
0,2,0,0;
0,0,2,0;
0,0,0,3];

stiff = 800000*[ 1,-1, 0, 0;
-1, 3,-2, 0;
0,-2, 5,-3;
0, 0,-3, 7];

damping = alpha1*mass + alpha2*stiff;

% Define earthquake record ...

Elcentro = [
0.00;
14.56; 13.77; 6.13; 3.73; 1.32; -6.81; -16.22; -22.41;
-26.67; -24.10; -21.86; -14.02; -1.62; 1.67; 18.12; 39.26;
28.40; 16.72; 16.06; 12.54; 0.15; -8.60; -14.00; -12.01;
-3.34; 5.14; 14.04; 16.36; 29.90; 43.09; 32.13; 17.28;
18.44; 22.31; 13.24; -0.91; -6.81; -15.20; -14.88; -9.54;
-13.32; -10.60; 11.80; 10.20; -12.90; -22.88; 1.42; 21.73;
7.92; 6.47; 11.39; -8.29; -33.44; -26.25; -23.99; -31.05;
-25.42; -26.97; -38.48; -35.01; -27.78; -33.35; -29.56; -13.60;
-13.97; -11.63; -1.31; 13.86; 35.71; 27.62; 25.33; 46.45;
44.12; 28.02; 19.13; 12.57; 11.60; 12.29; 9.73; 10.38;
18.93; 29.43; 30.25; 26.23; 23.15; 22.29; 21.11; 14.84;
11.02; 14.74; 21.65; 26.66; 28.45; 30.34; 30.21; 28.82;
31.31; 26.66; 15.29; 9.10; 3.58; -6.65; -15.33; -18.05;
-17.46; -11.68; -3.48; -5.25; -15.43; -25.38; -32.32; -29.67;
-22.01; -17.97; -22.34; -30.51; -38.89; -45.29; -55.26; -67.84;
-76.74; -84.58; -86.83; -86.63; -80.87; -69.83; -59.81; -54.51;
-48.94; -39.94; -38.67; -40.13; -38.04; -34.78; -31.30; -21.29;
-8.85; 3.58; 13.97; 30.88; 35.09; 8.84; -10.41; -16.30;
-20.10; -3.69; 33.71; 51.69; 43.52; 27.04; 34.70; 38.30;
36.60; 29.85; -9.22; -23.98; -17.28; -14.20; 11.58; 47.99;
64.26; 70.36; 68.93; 62.67; 45.30; 28.43; 23.73; 28.51;
34.46; 43.12; 43.39; 26.64; 21.23; 31.23; 40.30; 16.76;
-10.80; -22.88; -0.51; 12.51; -10.01; -12.52; -6.05; -24.18;
-43.17; 7.01; 50.06; 54.19; 68.83; 69.26; 57.53; 8.73;
-38.98; -43.59; -51.12; -41.86; -19.51; -15.78; -11.54; -2.18;
-2.19; 20.84; 23.40; -0.91; -11.45; -22.94; -40.84; -53.17;
-44.02; -21.96; -3.06; -0.49; -9.59; -11.83; -9.86; -7.02;
0.88; 24.00; 50.57; 56.01; 51.75; 52.30; 47.79; 16.77;
-6.09; -20.13; -20.42; -4.35; -0.33; -6.33; 1.74; 19.90;
37.97; 53.36; 57.68; 53.32; 28.37; 7.54; 6.90; 4.73;
-4.97; -21.97; -28.16; -23.25; -28.43; -18.57; -3.93; -16.60;
-31.92; -35.08; -34.77; -30.85; -30.37; -6.83; 24.57; 27.39;
23.72; 29.30; 5.79; -35.35; -54.31; -57.00; -59.38; -19.00;
32.14; 7.01; -21.44; -21.26; -26.01; -35.66; -33.68; -24.09;
-21.39; -26.29; -27.73; -23.91; -14.27; -6.41; -0.59; 4.21;
1.06; -4.59; -15.16; -18.12; -7.09; 2.41; 6.08; 1.20;
-10.53; -8.69; -4.70; -13.05; -9.51; -0.67; -7.10; -7.19;
5.36; 9.80; 7.87; 12.63; 20.28; 22.44; 24.86; 24.31;
15.25; 6.54; 6.86; 11.08; 18.12; 23.64; 22.47; 23.76;
21.57; 19.80; 25.47; 27.94; 14.81; -6.94; -13.57; -7.80;
-2.65; -3.97; -1.01; 0.20; 5.14; 2.31; -16.89; -18.18;
-11.55; -12.10; -3.07; 6.16; 8.87; 19.72; 21.31; 16.45;
12.70; 4.15; -1.59; 0.57; 11.04; 21.39; 27.42; 15.82;
-9.25; -23.67; -17.34; 0.67; 8.69; 11.01; 12.01; 11.12;
13.62; 11.31; 9.71; 17.68; 15.83; 5.59; 0.60; -2.53;
-4.46; -4.09; -7.81; -15.06; -19.43; -15.49; -8.93; -5.44;
-1.33; -1.35; 7.19; 12.52; 0.21; -4.89; 0.22; 11.51;
20.17; 18.62; 16.72; 13.74; 10.62; 11.14; 12.74; 13.34;
14.16; 14.59; 13.81; 10.94; 5.91; 5.74; 14.53; 19.56;
9.27; -6.75; -14.69; -14.25; -11.48; -3.96; 6.54; 10.59;
2.26; -12.97; -24.01; -25.73; -18.75; -6.99; -2.22; 0.17;
3.63; -1.70; -4.73; 6.13; 20.26; 22.36; 19.58; 9.14;
-0.09; -0.23; -1.00; 4.82; 15.39; 19.75; 7.45; -0.99;
0.77; -1.64; -9.34; -14.34; -13.86; -6.61; 5.49; 6.30;
2.08; -2.88; -2.25; 2.99; 3.30; -1.57; -6.86; -8.44;
-11.03; -16.68; -17.43; -12.39; -12.20; -15.43; -16.28; -17.99;
-20.21; -21.58; -23.05; -23.79; -23.64; -21.60; -15.32; -5.72;
11.14; 28.30; 33.22; 30.27; 25.49; 18.45; 7.62; -1.66;
-7.86; -10.04; -5.71; -5.16; -12.56; -19.39; -22.16; -18.18;
-12.00; -3.60; 3.99; 3.89; 0.85; 4.01; 6.05; 0.68;
-0.14; 1.50; 3.17; 4.96; -0.03; -7.20; -4.85; -0.60;
-3.05; -6.01; -6.08; -6.40; -6.70; -5.70; -6.14; -8.41;
-10.05; -12.35; -15.72; 0.00 ];

size(Elcentro);

% plot ground motion

plot( [ 0:0.02:10 ] , Elcentro)
axis([ 0, 12, -100, 100]);

max(abs(Elcentro));

% setup arrays for displacement, velocity, acceleration ....

velocity = zeros(4,1);
displacement = zeros(4,1);
acceleration = zeros(4,1);

% Create response array .....

response = zeros( nsteps + 1, 2);

% Compute M-hat matrix and r matrix .....

r = [ 1; 1; 1; 1 ];
mhat = mass+0.5*.02*damping+.25*(.02)^2*stiff;

% Begin main loop for Newmark iteration ....

for i = 1:nsteps

fprintf('Starting Loop %3.0f\n', i);

% P vector ....

if i <= 500
accel = Elcentro(i,1);
else
accel = 0;
end

P = -mass*r*accel;

% Compute P-hat array (see eqn 4.10).

phat = P-damping*velocity-stiff*displacement-dt*stiff*velocity-dt*((1-.5)*damping+(dt/2)*(1-2*.25)*stiff)*acceleration;

% Compute acceleration at new timestep (see eqn 4.8)

new_acceleration = mhat\phat;

% Compute new velocity (eqn 4.6)

new_velocity = velocity+dt*((1-.5)*acceleration+.5*new_acceleration);

% Compute new displacement (eqn 4.7)

new_displacement = displacement+dt*velocity+(dt^2/2)*((1-2*.25)*acceleration+2*.25*new_acceleration);

% Update velocity, acceleration, displacement .....

velocity = new_velocity;
displacement = new_displacement;

% Store response .....

response(i+1,1) = i*dt;
response(i+1,2) = displacement(1,1);
response(i+1,3) = velocity(1,1);
response(i+1,4) = acceleration(1,1);
end

% plot response

plot(response (:,1), response(:,2))
title('Displacement vs. Time for Top Floor')
xlabel('Time(s)')
ylabel('Displacement(m)')
plot(response (:,1), response(:,3))
title('Velocity vs. Time for Top Floor')
xlabel('Time(s)')
ylabel('Velocity(m/s)')
plot(response (:,1), response(:,4))
title('Acceleration vs. Time for Top Floor')
xlabel('Time(s)')
ylabel('Acceleration(m/s/s)')

%find max displacment
maxdisp=max(abs( response(:,2)))

%print shape of structure

if time==2|time==4|time==6|time==8|time==10|time==12,
new_displacement;

end

%if building will collapse

%if max disp/ height <.015 then building will not collapse

collapse_parameter=(maxdisp)/(120)