clear all
clc
data=dir('*.txt')
N=length(data);
for(k=1:N)
filename = data(k).name;% an acceleration time history
% the following block of code reads the specified file and integrates it to
% determine the velocity time history
file_in = fopen(filename, 'r');
% [record_dt ans] = fscanf(file_in, '%f', 1);
% ans = fgetl(file_in);
[record, np] = fscanf(file_in, '%f');
fclose(file_in);
%%%%%
id = findstr('.',filename);
sname=[filename(1:id(1)-4)];
% id1 = findstr(' ',filename);
% site=[filename(id1(2)+1:id1(2)+2)];
%%%%%
% Dt=record(1);
% Accelerate=record(2:end);
Dt=0.005; %%%可以修改Dt的值
Accelerate=record;
count=length(Accelerate);
time=0:Dt:(count-1)*Dt; %单位 s
% ***********精确法计算各反应***********
%初始化各储存向量
Displace=zeros(1,count); %相对位移
Velocity=zeros(1,count); %相对速度
AbsAcce=zeros(1,count); %绝对加速度
% ***********A,B矩阵***********
DampA=[0.05]; %三个阻尼比
% TA=[0.01 0.02 0.03 0.05 0.075 0.1 0.15 0.2 0.25 0.3 0.4 0.5 0.75 1 1.5 2 3 4 5 6 7.5 10]; %TA=0.000001:0.02:6; %结构周期
TA=[3.0]; % 只计算一个周期点
%记录计算得到的反应,MDis为某阻尼时最大相对位移,MVel为某阻尼
%时最大相对速度,MAcc某阻尼时最大绝对加速度,用于画图
MDis=zeros(3,length(TA));
MVel=zeros(3,length(TA));
MAcc=zeros(3,length(TA));
% figure
j=1; %在下一个循环中控制不同的阻尼比
for Damp=[0.05]
t=1; %在下一个循环中控制不同的结构自振周期
for T=TA
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 i=1:(count-1) %根据地震记录,计算不同的反应
Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i+1);
Velocity(i+1)=A(2,1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i+1);
AbsAcce(i+1)=-2*Damp*Frcy*Velocity(i+1)-Frcy^2*Displace(i+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
if mod(k,2)
xa=AbsAcce;
xv=Velocity;
xd=Displace;
else
ya=AbsAcce;
yv=Velocity;
yd=Displace;
figure
plot(xa,ya,'k-','linewidth',2)
hold on
xlabel('PSA x direction (cm/s^2)','FontSize',14,'FontWeight','Demi');
ylabel('PSA y direction (cm/s^2)','FontSize',14,'FontWeight','Demi');
set(gca,'fontsize',12,'fontweight','demi');
set(gca,'linewidth',2);
set(gca,'ticklength',[0.02;0.05]);
axis([-100 100 -100 100])
set(gca,'xtick',[-80:20:80]);
set(gca,'ytick',[-80:20:80]);
set(gca,'xticklabel',[-80:20:80]);
set(gca,'yticklabel',[-80:20:80]);
% axis([-50 50 -50 50])
for angle=1:180
rt=xa*cos(angle/180*pi)+ya*sin(angle/180*pi);
rth(angle)=max(abs(rt));
rx(angle)=rth(angle)*cos(angle/180*pi);
ry(angle)=rth(angle)*sin(angle/180*pi);
rx1(angle)=rth(angle)*cos(angle/180*pi+pi);
ry1(angle)=rth(angle)*sin(angle/180*pi+pi);
end
plot(rx,ry,'k-','linewidth',3)
plot(rx1,ry1,'k--','linewidth',3)
RotD_all(:,k/2)=rth';
RotD100(k/2,1)=max(rth);
RotD00(k/2,1)=min(rth);
RotD50(k/2,1)=median(rth);
id=find(rth==max(rth));
Rot_angle(k/2,1)=id;
x1=[-max(rth)*cos(id/180*pi) max(rth)*cos(id/180*pi)];
y1=[-max(rth)*sin(id/180*pi) max(rth)*sin(id/180*pi)];
plot(x1,y1,'r-','linewidth',2.5)
if ~exist('SA30')
mkdir('SA30')
end
paths=[pwd,'\SA30\']; % cd返回当前目录(pwd)
saveas(gcf,[paths,sname,'_ra.jpg']);
saveas(gcf,[paths,sname,'_ra.fig']);
close
figure
plot(xv,yv,'k-','linewidth',2)
xlabel('Response velocity x direction (cm/s)','FontSize',14,'FontWeight','Demi');
ylabel('Response velocity y direction (cm/s)','FontSize',14,'FontWeight','Demi');
set(gca,'fontsize',12,'fontweight','demi');
set(gca,'linewidth',2);
set(gca,'ticklength',[0.02;0.05]);
axis([-100 100 -100 100])
set(gca,'xtick',[-100:20:100]);
set(gca,'ytick',[-100:20:100]);
set(gca,'xticklabel',[-100:20:100]);
set(gca,'yticklabel',[-100:20:100]);
saveas(gcf,[paths,sname,'_rv.jpg']);
saveas(gcf,[paths,sname,'_rv.fig']);
close
end
% Displace=zeros(1,count);%初始化各储存向量,避免下次不同周期计算时引用到前一个周期的结果
% Velocity=zeros(1,count);
% AbsAcce=zeros(1,count);
% t=t+1;
end
end
% fid=fopen([paths,strcat(sname,'.txt')],'w');
% fprintf(fid,'%5.2f %5.2f %5.2f \n',MAcc);
% fclose(fid);
%
% ra5(k,1)=MAcc(1,8);%T=0.2周期谱加速度值
% ra5(k,2)=MAcc(1,12);%T=0.5周期谱加速度值
% ra5(k,3)=MAcc(1,14);%T=1.0周期谱加速度值
% ra5(k,4)=MAcc(1,17);%T=3.0周期谱加速度值
str{k}=sname;
end;
% close all
% save Resp_acc.txt ra5 -ascii
fid = fopen('Name.dat', 'wt');
for i=1:N/2
fprintf(fid, '%s\n', str{i*2});
end