%%%% the program is to compute the parameter of PEER Ground Motion
clear
clc
close all
filename=dir('over50\*.AT2');
n=length(filename);
%========= the Names of items which we want to get ============%
File_Name=cell(n,1);PGA=zeros(n,1);PGV=zeros(n,1);PGD=zeros(n,1);
Uniform_Dur=zeros(n,1);Bracketed_Dur=zeros(n,1);D90=zeros(n,1);
AI=zeros(n,1);MIV=zeros(n,1);MID=zeros(n,1);CAV=zeros(n,1);CAVSTD=zeros(n,1);
SA_02=zeros(n,1);SA_1=zeros(n,1);SV_02=zeros(n,1);SV_1=zeros(n,1);SD_02=zeros(n,1);
SD_1=zeros(n,1);EPA=zeros(n,1);EPV=zeros(n,1);VSI=zeros(n,1);
g=980;
%=======================================================================%
for k=1:n
%=========== read the namelist of files =======================%
ss=length(filename(k).name);
File_Name{k}=filename(k).name(1:ss-4);
% to read the kth file
[Infor,Acceleration]=Read_Files_Peer(['over50\',filename(k).name]);
Dt=Infor.Dt;
Acceleration=Acceleration*g; % the unit of acceleration is gal
count=length(Acceleration);
Velocity=zeros(count,1);
Displacement=zeros(count,1);
%============ COMPUTE PGA , PGV and PGD ========================%
PGA(k)=max(abs(Acceleration));
for ii=1:(count-1)
Velocity(ii+1)=Velocity(ii)+(Acceleration(ii)+Acceleration(ii+1))*Dt*(1/2);
Displacement(ii+1)=Displacement(ii)+Velocity(ii)*Dt+(Acceleration(ii)*(1/3)+Acceleration(ii+1)*(1/6))*Dt*Dt;
end
PGV(k)=max(abs(Velocity));
PGD(k)=max(abs(Displacement));
%============== COMPUTE Uniform Duration : Uniform_Dur ===============%
Default1=0.05*PGA(k);
n1=0;
for k1=1:count
if abs(Acceleration(k1))>Default1
n1=n1+1;
end
end
Uniform_Dur(k)=n1*Dt;
%=========== COMPUTE Bracketed Duration : Bracketed_Dur ===============%
Default2=Default1;
C=find(abs(Acceleration)>Default2);
BT1=max(C);
BT2=min(C);
Bracketed_Dur(k)=(BT1-BT2+1)*Dt;
%=============== COMPUTE Significant Duration : D90=============%
IT=0;
for k2=1:count-1
IT=IT+(Acceleration(k2)*Acceleration(k2)+Acceleration(k2+1)*Acceleration(k2+1))*Dt*(1/2);
end
ITt=0;
for k3=1:count-1
ITt=ITt+(Acceleration(k3)*Acceleration(k3)+Acceleration(k3+1)*Acceleration(k3+1))*Dt*(1/2);
if ITt>0.05*IT
t1=(k3-1)*Dt;
break;
end
end
for k4=1:count-1
ITt=ITt+(Acceleration(k4)*Acceleration(k4)+Acceleration(k4+1)*Acceleration(k4+1))*Dt*(1/2);
if ITt>0.95*IT
t2=(k4-1)*Dt;
break;
end
end
D90(k)=t2-t1;
%============ COMPUTE Arias Intensity : AI =======================%
AI(k)=IT*pi/(2*g);
%= COMPUTE Maximum Increment Velocity and Displacement : MIV and MID =%
IV=zeros(count,1);ID=zeros(count,1);
j=1;jj=1;k5=1;
while j+k5<=count
if Acceleration(j)*Acceleration(j+k5)>0
k5=k5+1;
else
IV(jj)=sum(Acceleration(j:j+k5-1))*Dt;
jj=jj+1;
j=j+k5;
k5=1;
end
end
j=1;jj=1;k6=1;
while j+k6<=count
if Acceleration(j)*Acceleration(j+k6)>0
k6=k6+1;
else
ID(jj)=sum(Velocity(j:j+k6-1))*Dt;
jj=jj+1;
j=j+k6;
k6=1;
end
end
MIV(k)=max(abs(IV));%unit:cm/s
MID(k)=max(abs(ID));%unit:cm
%============= COMPUTE Cumulative Absolute Velocity : CAV ===========%
CAV_t=0;
ABS_ACC=abs(Acceleration);
for k7=1:count-1
CAV_t=CAV_t+(ABS_ACC(k7)+ABS_ACC(k7+1))*Dt*(1/2);
end
CAV(k)=CAV_t;
%====== COMPUTE Standard Cumulative Absolute Velocity : CAVSTD =======%
CAVSTD_t=0;
a=count*Dt;
if fix(a)==a
a=count*Dt;
else
a=fix(a)+1;
Acceleration((count+1):fix(a/Dt))=0;
end
for t=1:a
for k3=fix((t-1)/Dt)+1:1:fix(t/Dt)
if max(abs(Acceleration((fix((t-1)/Dt)+1):1:fix(t/Dt))))>=24.5
CAVSTD_t=CAVSTD_t+abs(Acceleration(k3))*Dt;
continue;
else
break;
end
end
end
CAVSTD(k)=CAVSTD_t;
%==COMPUTE Response Spectra at 0.2 sec : SA_02 , SV_02 , SD_02 ======%
%==COMPUTE Response Spectrua at 1 sec : SA_1 , SV_1 , SD_1 ==========%
%==COMPUTE Effective Peak Acceleration and Velocity : EPA ,EPV ======%
time=0:Dt:Dt*(count-1);
Displace=zeros(count,1);Veloci=zeros(count,1);AbsAcce=zeros(count,1);
TA=[0 0.01 0.02 0.025 0.03 0.032 0.034 0.036 0.038 0.04 0.042 0.044 ...
0.046 0.048 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.085 0.09...
0.095 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.22...
0.24 0.25 0.26 0.28 0.3 0.32 0.34 0.35 0.36 0.38 0.4 0.42 0.44...
0.45 0.46 0.48 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.1...
1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.2 2.4 2.5 2.6 2.8 3 3.2 3.4 3.6...
3.8 4 4.2 4.4 4.6 4.8 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10];
MDis=zeros(length(TA),1);
MVel=zeros(length(TA),1);
MAcc=zeros(length(TA),1);
t=1;
for T=[0 0.01 0.02 0.025 0.03 0.032 0.034 0.036 0.038 0.04 0.042 0.044 ...
0.046 0.048 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.085 0.09...
0.095 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.22...
0.24 0.25 0.26 0.28 0.3 0.32 0.34 0.35 0.36 0.38 0.4 0.42 0.44...
0.45 0.46 0.48 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.1...
1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.2 2.4 2.5 2.6 2.8 3 3.2 3.4 3.6...
3.8 4 4.2 4.4 4.6 4.8 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10]
Fracy=2*pi/T;
Damp=0.05;
dt=Dt;% just for convenience
DamFracy=Fracy*sqrt(1-Damp*Damp);
e_t=exp(-Damp*Fracy*Dt);
s=sin(DamFracy*Dt);
c=cos(DamFracy*Dt);
A=zeros(2,2);
A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);
A(1,2)=e_t*s/DamFracy;
A(2,1)=-Fracy*e_t*s/sqrt(1-Damp*Damp);
A(2,2)=e_t*(c-Damp*s/sqrt(1-Damp*Damp));
d2t=(2*Damp*Damp-1)/(Fracy*Fracy*Dt);
d3t=Damp/(Fracy*Fracy*Fracy*Dt);
B=zeros(2,2);
B(1,1)=e_t*((d2t+Damp/Fracy)*s/DamFracy+(2*d3t+1/Fracy^2)*c)-2*d3t;
B(1,2)=-e_t*(d2t*s/DamFracy+2*d3t*c)-1/Fracy^2+2*d3t;
B(2,1)=e_t*((d2t+Damp/Fracy)*(c-Damp*s/sqrt(1-Damp*Damp))-(2*d3t+1/Fracy^2)*(DamFracy*s+Damp*Fracy*c))+1/(Fracy^2*Dt);
B(2,2)=e_t*(1/(Fracy^2*Dt)*c+s*Damp/(Fracy*DamFracy*Dt))-1/(Fracy^2*Dt);
for ii=1:count-1
Displace(ii+1)=A(1,1)*Displace(ii)+A(1,2)*Veloci(ii)+B(1,1)*Acceleration(ii)+B(1,2)*Acceleration(ii+1);
Veloci(ii+1)=A(2,1)*Displace(ii)+A(2,2)*Veloci(ii)+B(2,1)*Acceleration(ii)+B(2,2)*Acceleration(ii+1);
AbsAcce(ii+1)=-2*Damp*Fracy*Veloci(ii+1)-Fracy^2*Displace(ii+1);
end
MDis(t)=max(abs(Displace));
MVel(t)=max(abs(Veloci));
if T==0
MAcc(t)=max(abs(Acceleration));
else
MAcc(t)=max(abs(AbsAcce));
end
t=t+1;
Displace=zeros(count,1);
Veloci=zeros(count,1);
AbsAcce=zeros(count,1);
end
SA_02(k)=MAcc(35);%====== 0.2s
SA_1(k)=MAcc(63);%======= 1.0s
SV_02(k)=MVel(35);
SV_1(k)=MVel(63);
SD_02(k)=MDis(35);
SD_1(k)=MDis(63);
EPA(k)=mean(MAcc(25:53))/2.5;% the interval is 0.1-0.5s
EPV(k)=mean(MVel(59:65))/2.5;% the interval is 0.8-1.2s
%=========COMPUTE Housner Velocity Spectrum Intensity : VSI =======%
VSI_t=0;
for kk=25:76
VSI_t=VSI_t+(MVel(kk)+MVel(kk+1))*(TA(kk+1)-TA(kk))/2;
end
VSI(k)=VSI_t;
end
Columns={'File_Name','PGA','PGV','PGD','Uniform_Dur','Bracketed_Dur','