地震波转谱程序

fid = fopen（‘D:\世茂胶南\安评地震\画图用\S645-5.txt‘）;
[Acceleratecount] = fscanf（fid‘%g‘）; %count 读入的记录的量
Accelerate=Accelerate*40/1616; %单位统一为 m和s
time=0:0.02:（count-1）*0.02; %单位 s
Displace=zeros（1count）; %相对位移
Velocity=zeros（1count）; %相对速度
AbsAcce=zeros（1count）; %绝对加速度
% ***********AB矩阵***********
DampA=[00.050.1]; %三个阻尼比
TA=0.0:0.05:6; %TA=0.000001:0.02:6; %结构周期
Dt=0.02; %地震记录的步长
%记录计算得到的反应，MDis为某阻尼时最大相对位移，MVel为某阻尼
%时最大相对速度，MAcc某阻尼时最大绝对加速度，用于画图
MDis=zeros（3length（TA））;
MVel=zeros（3length（TA））;
MAcc=zeros（3length（TA））;

j=1; %在下一个循环中控制不同的阻尼比
for Damp=[00.050.1]

t=1; %在下一个循环中控制不同的结构自振周期
for T=0.0:0.05: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（22）;
A（11）=e_t*（s*Damp/sqrt（1-Damp*Damp）+c）;
A（12）=e_t*s/DamFrcy;
A（21）=-Frcy*e_t*s/sqrt（1-Damp*Damp）;
A（22）=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（22）;
B（11）=e_t*（（d_f+Damp/Frcy）*s/DamFrcy+（2*d_3t+1/Frcy^2）*c）-2*d_3t;
B（12）=-e_t*（d_f*s/DamFrcy+2*d_3t*c）-1/Frcy^2+2*d_3t;
B（21）=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（22）=e_t*（1/（Frcy^2*Dt）*c+s*Damp/（Frcy*DamFrcy*Dt））-1/（Frcy^2*Dt）;

for i=1:（count-1） %根据地震记录计算不同的反应
Displace（i+1）=A（11）*Displace（i）+A（12）*Velocity（i）+B（11）*Accelerate（i）+B（12）*Accelerate（i+1）;
Velocity（i+1）=A（21）*Displace（i）+A（22）*Velocity（i）+B（21）*Accelerate（i）+B（22）*Accelerate（i+1）;
AbsAcce（