! 平均加速度法计算单自由度SDOF系统地震响应计算源程序
! Sourse code for seismic response of SDOF system with average acceleration method.
! 输入数组EQ(:,2) 地震加速度记录，时间，加速度
! Input Array: EQ(:,2) Earthquake acceleration record: Time, Acceleration
! 输出数组 Disp(:)
! Output Array: Disp(:)
! 程序编写：陆新征，清华大学土木工程系，luxinzheng@sina.com
! By Xinzheng LU, Civil Engineering of Tsinghua University, luxinzheng@sina.com

program Main
real*8 EQ(4000,2),Disp(4000)
real*8 Te, Damp

Te=0.78219 ! Period of system

Damp=0.05d0; ! Damping ratio

open(55,file="EQ01.txt")
do I=1, size(EQ,1)
read(55,*) EQ(I,:)
end do
close(55)

call Average_Acce(EQ,Disp,Te, Damp,size(EQ,1))

open(66,file='disp.txt')
do I=1, size(EQ,1)
write(66,*) EQ(I,1), Disp(I)
end do
close(55)

stop
end program

subroutine Average_Acce(EQ,Disp,Te, Damp,N)
implicit none
integer :: N
real*8 EQ(N,2)
real*8 Disp(N)
real*8 Te, Damp
real*8 k, m, c, Omega,pi
real*8 a1,a2,v1,v2,d1,d2
real*8 dt
integer I,J, STEP

STEP=100 ! Number of substep for integration

pi=4.d0*atan(1.d0)

m=1.0;
K=4.d0*pi*pi*m/Te/Te;

Omega=2.d0*pi/Te

c=2.d0*Damp*sqrt(k*m);

a1=0.d0; a2=0.d0;
v1=0.d0; v2=0.d0;
d1=0.d0; d2=0.d0;
Disp(1)=0.d0;

do I=2, N
dt=(EQ(I,1)-EQ(I-1,1))/real(STEP);
do J=1, STEP
a1=a2; v1=v2; d1=d2;

a2=EQ(I-1,2)+(EQ(I,2)-EQ(I-1,2))*real(J)/real(STEP)
a2=a2+c/m*(v1+0.5d0*a1*dt)+k/m*(d1+v1*dt+0.25d0*a1*dt*dt);
a2=-a2/(1.d0+0.5d0*c/m*dt+0.25d0*k/m*dt*dt)

v2=v1+0.5d0*(a1+a2)*dt;
d2=d1+v1*dt+0.25d0*(a1+a2)*dt*dt
end do
Disp(I)=d2
end do

return
end subroutine