clear
clc
close all
filename=dir('example_liang\*.AT2');
n=length(filename);
for k=1:n
fid=fopen(['example_liang\',filename(k).name],'r');
%%%%%%%%%%%%%%%to find the time step Dt
text=textscan(fid,'%s',4,'delimiter','\n');% read the data at the fist four lines
a=text{1,1}{4,1};
b=a(21:25);
Dt=str2double(b);
%%%%%%%%%%%%%%to read the strong motion data
Accelerate=textread(['example_liang\',filename(k).name],'%f','headerlines',4);
count=length(Accelerate);
time=0:Dt:(count-1)*Dt;% I am not sure the use of Dt
Displace=zeros(1,count);
Velocity=zeros(1,count);
AbsAcce=zeros(1,count);
DampA=[0,0.05,0.1];%%%The Damping facor can be modified
TA=0:Dt:6;
MDis=zeros(3,length(TA));
MVel=zeros(3,length(TA));
MAcc=zeros(3,length(TA));
j=1;
for Damp=[0.005,0.05,0.1]
t=1;
for T=0:Dt:6
Frcy=2*pi/T;
DamFrcy=Frcy*sqrt(1-Damp*Damp);
e_t=exp(-Damp*Frcy*Dt);
s=sin(DamFrcy*Dt);
c=cos(DamFrcy*Dt);
A=zeros(2,2);
A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);
A(1,2)=e_t*s/DamFrcy;
A(2,1)=-Frcy*e_t*s/sqrt(1-Damp*Damp);
A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);
d_f=(2*Damp^2-1)/(Frcy^2*Dt);
d_3t=Damp/(Frcy^3*Dt);
B=zeros(2,2);
B(1,1)=e_t*((d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy^2)*c)-2*d_3t;
B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy^2+2*d_3t;
B(2,1)=e_t*((d_f+Damp/Frcy)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/Frcy^2)*(DamFrcy*s+Damp*Frcy*c))+1/(Frcy^2*Dt);
B(2,2)=e_t*(1/(Frcy^2*Dt)*c+s*Damp/(Frcy*DamFrcy*Dt))-1/(Frcy^2*Dt);
for ii=1:(count-1)
Displace(ii+1)=A(1,1)*Displace(ii)+A(1,2)*Velocity(ii)+B(1,1)*Accelerate(ii)+B(1,2)*Accelerate(ii+1);
Velocity(ii+1)=A(2,1)*Displace(ii)+A(2,2)*Velocity(ii)+B(2,1)*Accelerate(ii)+B(2,2)*Accelerate(ii+1);
AbsAcce(ii+1)=-2*Damp*Frcy*Velocity(ii+1)-Frcy^2*Displace(ii+1);
end
MDis(j,t)=max(abs(Displace));
MVel(j,t)=max(abs(Velocity));
if T==0.0
MAcc(j,t)=max(abs(Accelerate));
else
MAcc(j,t)=max(abs(AbsAcce));
end
Displace=zeros(1,count);
Velocity=zeros(1,count);
AbsAcce=zeros(1,count);
t=t+1;
end
j=j+1;
end
h=figure(k);
plot(TA,MAcc(1,:),'-.b',TA,MAcc(2,:),'-r',TA,MAcc(3,:),':k');
xlabel('Tn(s)');
ylabel('absolute acceleration(g)');
legend('汎=0','汎=0.05','汎=0.1');
grid
ss=length(filename(k).name);
subs=filename(k).name(1:ss-4);
plus_Acceleration='_Acceleration';
subs=[subs,plus_Acceleration];
title(subs);
set(gca,'xscale','log');
saveas(h,[subs,'.fig']);
close(h);
end