%ce445_Response_Spec_anim.m
%Animation function for the lateral vibration of an SDOF building
% for  demonstration of response spectra
clear all; close all; clc; 
warning off all;

anim=0;
fs_lab=14; % font size for figure labels
%fs_lab1=14;  % font size for spectrum 

dind=0; % demo index, for subplotting
specfigno=1; % figure number

%%  GET USER INPUT and READ DATA
name='Input Parameters';
p1='damping ratio ';
p2='shortest Tn (sec)';
p3='longest Tn (sec)';
prompt={p1,p2,p3}; %,p4,p5};
numlines=1;
defaultanswer={'0.05','0.01', '10'};
answer=inputdlg(prompt, name, numlines,defaultanswer);
zeta=str2double(answer{1}); 
Tn1=str2double(answer{2});   
Tn2=str2double(answer{3});   
Tnvec=logspace(log10(Tn1),log10(Tn2),100);
Tnvec_demoSpec=[0.5,1:10];
Tnvec=unique([Tnvec,Tnvec_demoSpec]);
ndemoSpec=length(Tnvec_demoSpec);

colorvec=['r' ,'g','b','k','m','y','c'];
colorvec=[colorvec,colorvec]; %['r' ,'g','b','k','m','y','c'];
colorvec=[colorvec,colorvec]; 

load data     % dt and acc(in/sec^2)
N=length(acc);         
t=(0:N-1)*dt;       % time vector
f=386*acc;  % in/sec^2
wnvec=2*pi./Tnvec;
nf=length(wnvec);

%% STATE_SPACE MODEL
B = [ 0 ; -1 ]; C = [ 1  0 ]; D = 0; 
xo=0;
vo=0;
init=[xo,vo];

Sdvec=zeros(size(Tnvec));
for j=1:nf;
wn=wnvec(j);
Tn=2*pi/wn; % current Tn 
A = [ 0  1 ; -wn^2  -2*zeta*wn ]; % current A
sys=ss(A, B, C, D);

% compute response
x1=lsim(sys,f,t,init); 
Sd1=max(abs(x1)); % peak displacement
tloc=find(abs(x1)==Sd1);

Sdvec(j)=Sd1;

if ismember(Tn, Tnvec_demoSpec)
dind=dind+1; % demo ind
col=colorvec(dind); % color

figure(specfigno)
subplot 323
hold off
plot(t,x1)
hold on;
plot(t(tloc),x1(tloc),'o','markerFaceColor',col); % put a circle
xlabel('t (s)','fontsize',fs_lab),
ylabel('x (in)','fontsize',fs_lab); 
ylim(1.5*Sd1*[-1 1]);
str=['T_n=',num2str(Tn),' s']; %, S_d= ', num2str(Sd1,4),' in'];
title(str,'fontsize',14,'backgroundcolor','y')
str=['S_d= ', num2str(Sd1,4),' in'];
text(t(tloc)+2,1.2*x1(tloc),str,'color','r');%,'backgroundcolor','y')


hold on
subplot(322)
hold on
plot(Tn,Sd1,'o','markerFaceColor',col); % 
ylabel('S_d(in)','fontsize',fs_lab); 
xlim([0.0001 Tnvec(end)])

subplot(324)
hold on
plot(Tn,2*pi/Tn*Sd1,'o','markerFaceColor',col); % 
ylabel('S_v(in/s)','fontsize',fs_lab); 
xlim([0.0001 Tnvec(end)])


subplot(326)
hold on;
plot(Tn,(2*pi/Tn)^2*Sd1,'o','markerFaceColor',col); % 
xlabel('Tn (s)','fontsize',fs_lab);
ylabel('S_a(in/s^2)','fontsize',fs_lab); 
xlim([0.0001 Tnvec(end)])
% title('S_a= \omega_n^2*S_d')


set(gcf, 'position',[100 100 800 800])
shg, 
pause(1)



end



end

% full spectrum

Svvec=(2*pi./Tnvec).*Sdvec;
Savec=(2*pi./Tnvec).^2.*Sdvec;


figure(specfigno);
subplot 322
plot(Tnvec,Sdvec);
str=['Response Spectra for El Centro Ground Motion, \zeta= ',num2str(zeta)];
text(-10,1.25*max(Sdvec),str,'fontsize',14)

subplot 324
plot(Tnvec,Svvec);

subplot 326
plot(Tnvec,Savec);

