>> fid=fopen('E:\工力所研究生\地震工程\作业\zuoye1\63.AT2');
>> ag=fscanf(fid,'%e');
>> Accelerate=ag(:,1);
maxA=max(abs(Accelerate))
Accelerate=0.5/maxA*Accelerate
>> N=length(Accelerate);
>> time=0:0.01:(N-1)*0.01;
>> Displace=zeros(1,N);
>> Velocity=zeros(1,N);
>> AbsAcce=zeros(1,N);
>>
>> Damp=0.05;
>> TA=0.0:0.01:6;
>> Dt=0.01;
>>
>> MaxD=zeros(1,length(TA));
>> MaxV=zeros(1,length(TA));
>> MaxA=zeros(1,length(TA));
>> t=1;
>> for T=0.0:0.01:6
NatualFrequency=2*pi/T;
DampFrequency=NatualFrequency*sqrt(1-Damp*Damp);
e_t=exp(-Damp*NatualFrequency*Dt);
s=sin(DampFrequency*Dt);
c=cos(DampFrequency*Dt);
A=zeros(2,2);
A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);
A(1,2)=e_t*s/DampFrequency;
A(2,1)=-NatualFrequency*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)/(NatualFrequency^2*Dt);
d_3t=Damp/(NatualFrequency^3*Dt);
B=zeros(2,2);
B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3t;
B(1,2)=e_t*(-d_f*s/DampFrequency-2*d_3t*c)-1/NatualFrequency^2+2*d_3t;
B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2)*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);
B(2,1)=-e_t*(((Damp/(NatualFrequency*Dt)+1)*s/DampFrequency)+(1/(NatualFrequency^2*Dt))*c)+1/(NatualFrequency^2*Dt);
B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualFrequency^2*Dt);
for i=1:(N-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*NatualFrequency*Velocity(i+1)-NatualFrequency^2*Displace(i+1);
end
MaxD(1,t)=max(abs(Displace));
MaxV(1,t)=max(abs(Velocity));
if T==0.0
MaxA(1,t)=max(abs(Accelerate));
else
MaxA(1,t)=max(abs(AbsAcce));
end
Displace=zeros(1,N);
Velocity=zeros(1,N);
AbsAcce=zeros(1,N);
t=t+1;
end
plot(TA,MaxA(1,:),'-b',TA,acc(:,1),'-r'),title('Acceleration'),xlabel('Tn(s)'),ylabel('Acceleration(g)'),legend('b=matlab','r=seismo'),grid;
plot(TA,980*MaxV(1,:),'-b',TA,vel(:,1),'-r'),title('Velocity'),xlabel('Tn(s)'),ylabel('Velocity(m/s)'),legend('b=matlab','r=seismo'),grid;
plot(TA,980*MaxD(1,:),'-b',TA,dis(:,1),'-r'),title('Displacement'),xlabel('Tn(s)'),ylabel('Displacement(m)'),legend('b=matlab','r=seismo'),grid;
SMaxA(1,:)=MaxA(1,:)/0.5
v=cumtrapz(time,Accelerate);
SMaxV(1,:)=MaxV(1,:)/max(abs(v));
d=cumtrapz(time,v);
SMaxD(1,:)=MaxD(1,:)/max(abs(d));
plot(TA,SMaxA(1,:),'-b',TA,acc(:,1)/0.5,'-r'),title('Acceleration'),xlabel('Tn(s)'),ylabel('Acceleration'),legend('b=matlab','r=seismo'),grid;
plot(TA,SMaxD(1,:),'-b',TA,dis(:,1)/980/max(abs(d)),'-r'),title('Displacement'),xlabel('Tn(s)'),ylabel('Displacement'),legend('b=matlab','r=seismo'),grid;
plot(TA,SMaxV(1,:),'-b',TA,vel(:,1)/980/max(abs(v)),'-r'),title('Velocity'),xlabel('Tn(s)'),ylabel('Velocity'),legend('b=matlab','r=seismo'),grid;
