地震动参数计算matlab程序

%%%% the program is to compute the parameter of PEER Ground Motion
clear
clc
close all
filename=dir（‘D:\matlab_my_data\work_for_liangyan‘）;
n=length（filename）;
%========= the Names of items which we want to get ============%
File_Name=cell（n1）;PGA=zeros（n1）;PGV=zeros（n1）;PGD=zeros（n1）;
Uniform_Dur=zeros（n1）;Bracketed_Dur=zeros（n1）;D90=zeros（n1）;
AI=zeros（n1）;MIV=zeros（n1）;MID=zeros（n1）;CAV=zeros（n1）;CAVSTD=zeros（n1）;
SA_02=zeros（n1）;SA_1=zeros（n1）;SV_02=zeros（n1）;SV_1=zeros（n1）;SD_02=zeros（n1）;
SD_1=zeros（n1）;EPA=zeros（n1）;EPV=zeros（n1）;VSI=zeros（n1）;
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
    [InforAcceleration]=Read_Files_Peer（[‘over50\‘filename（k）.name]）;
    Dt=Infor.Dt;
    Acceleration=Acceleration*g; % the unit of acceleration is gal
    count=length（Acceleration）;
    Velocity=zeros（count1）;
    Displacement=zeros（count1）;
    %============ 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（count1）;ID=zeros（count1）;
    j=1;jj=1;k5=1;
    while j+k5<=count
        if Acceleration（j）*Acceleration（j+k5）>0
            k5=k5+1;