% Lateral vibration of an SDOF system
% free or forced vibration

clear all; close all; clc; 
%%  GET USER INPUT
name='Input Parameters';
p1='damping ratio ';
p2='natural period (sec)';
p3='initial displacement (in)';
p4='initial velocity (in/s)';

prompt={p1,p2,p3,p4};
numlines=1;
defaultanswer={'0.05','1','2','0'};
answer=inputdlg(prompt, name, numlines,defaultanswer);
zeta=str2double(answer{1}); %use str2num if there is text
Tn=str2double(answer{2});   
xo=str2double(answer{3});  
vo=str2double(answer{4});  

% time vector
t=0:0.01:10*Tn;

% state space model
wn=2*pi/Tn;
A = [ 0  1 ; -wn^2  -2*zeta*wn ]; % current A
B = [ 0 ; -1 ]; C = [ 1  0 ]; D = 0; %these don’t change 
sys=ss(A, B, C, D);

% compute response
init=[xo,vo];
x1=lsim(sys,0*t,t,init); % this is in inches
Sd1=max(abs(x1)); % peak displacement

%% set geometric parameters
%h=wall heights W=floor width th=floor thickness
h=2*8;      W=2*.2;   th=2*.5;
slabx=[0,0,W,W]; 
slaby=[-th,0,0,-th];%plotting coordinates for floors
basex=[-W/5,-W/5,W+W/5,W+W/5]; 
basey=[-th,th,th,-th];%plotting coordinates for floors
%wall vertical coordinate arrays
s1=linspace(0,h,100);

% initialize plots
fs_lab=14;
figure('Units','Normalized','Position',[0.05,0.05,0.85,0.5]);
%deformation vs t
Ha_plot1=subplot('Position',[.55 .2 .4 .6]);
axis([t(1) t(end) -1.2*max(x1) 1.2*max(x1)]),
grid on; 
xlabel('t (s)','fontsize',fs_lab),
ylabel('x (in)','fontsize',fs_lab); 
tit1=['System: \zeta=', num2str(zeta),' and T_n= ', num2str(Tn),' s'];
tit2=['Initial conditions: xo =', num2str(xo),' in and vo = ', num2str(vo),' in/s'];
titstr={tit1;tit2};
title(titstr,'fontsize',fs_lab)


Hl_plotx1=line(t(1),x1(1),'Color','b','linewidth',2); 

%animation
Ha_anim=subplot('Position',[.1,.2,.4,.8]);
axis off; 
axis([-W 2*W 0 2*h]),% grid on; axis off;
%draw 1st floor walls
v1=x1(1)*(-2*s1.^3/h^3+3*s1.^2/h^2);
Hl_Lwall1=line(v1,s1,'Color','k','Linewidth',3);
Hl_Rwall1=line(v1+W,s1,'Color','k','Linewidth',3);

Hp_floor1=patch(slabx+x1(1),slaby+h,'b');
Hp_base=patch(basex,basey,'r');


%draw initial plots and hold
drawnow 

%Run animation
scale=0.05/Sd1;
tic;
for n=1:length(t)
    %plots
    set(Hl_plotx1,'XData',t(1:n),'YData',x1(1:n));
    v1=scale*x1(n)*(-2*s1.^3/h^3+3*s1.^2/h^2); % wall shapes
    % animation
    set(Hl_Lwall1,'XData',v1);
    set(Hl_Rwall1,'XData',v1+W);
    set(Hp_floor1,'XData',slabx+scale*x1(n));
     while toc<t(n), end %time delay
    drawnow
end

