module defpi
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    integer,parameter::md=5
end module defpi

program f90_acc_rng
    !*****************************************************************
    !Acceleration Response (time domein)
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i,j
    real(8),allocatable::acc(:)       !Acceleration (input data)
    real(8),allocatable::racc(:)      !Response acceleration (output data)
    real(8),allocatable::rvel(:)      !Response velocity (output data)
    real(8),allocatable::rdis(:)      !Response displacement (output data)
    real(8)::hh                       !Damping factor
    real(8)::Tp                       !Period
    real(8)::dt                       !Time increment
    character::fnameR*50              !Input & output file name
    character::dummy*50,linebuff*300
    real(8)::am,ac,ak,f1,f2,fm1,fm2,fmm
    real(8)::acc1,vel1,dis1
    real(8)::x1,y1,x2,y2,ddt
    real(8)::kx1,kx2,kx3,kx4
    real(8)::ky1,ky2,ky3,ky4

    call getarg(1,dummy);read(dummy,*) hh
    call getarg(2,dummy);read(dummy,*) Tp
    call getarg(3,fnameR)  !Input file name (time history data of acceleration)

    open(11,file=fnameR,status='old')
    read(11,'(a)') linebuff
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    allocate(acc(1:ndata),racc(1:ndata),rvel(1:ndata),rdis(1:ndata))
    do i=1,ndata
        read(11,*) acc(i)
    end do
    close(11)

    ddt=dt/dble(md)
    am=1.0D0                 !Mass
    ak=4.0D0*pi*pi*am/Tp/Tp  !Stiffness
    ac=2.0D0*hh*sqrt(ak*am)  !Damping constant

    acc1=0.0D0
    vel1=0.0D0
    dis1=0.0D0
    do i=1,ndata
        f1=0.0D0
        if(2<=i)f1=-am*acc(i-1)
        f2=-am*acc(i)
        x1=dis1
        y1=vel1
        do j=1,md
            fm1=f1+(f2-f1)*dble(j-1)/dble(md)
            fm2=f1+(f2-f1)*dble(j)/dble(md)
            fmm=0.5D0*(fm1+fm2)
            call FUNC(am,ac,ak,fm1,x1,y1,kx1,ky1)
            call FUNC(am,ac,ak,fmm,x1+ddt*kx1/2.0D0,y1+ddt*ky1/2.0D0,kx2,ky2)
            call FUNC(am,ac,ak,fmm,x1+ddt*kx2/2.0D0,y1+ddt*ky2/2.0D0,kx3,ky3)
            call FUNC(am,ac,ak,fm2,x1+ddt*kx3,y1+ddt*ky3,kx4,ky4)
            x2=x1+ddt/6.0D0*(kx1+2.0D0*kx2+2.0D0*kx3+kx4)
            y2=y1+ddt/6.0D0*(ky1+2.0D0*ky2+2.0D0*ky3+ky4)
            x1=x2
            y1=y2
        end do
        acc1=f2/am-ac/am*y2-ak/am*x2
        vel1=y2
        dis1=x2
        racc(i)=acc1+acc(i)   !Absolute acceleration
        rvel(i)=vel1          !Relative velocity
        rdis(i)=dis1          !Relative displacement
    end do

    write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'dt,',dt; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'ndata,',ndata; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'Tp,',Tp; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(linebuff,*) 'dt,',dt; call del_spaces(linebuff);write(6,'(a)') trim(linebuff)
    write(6,'(a)') 'time(sec),acc(gal),racc(gal),rvel(kine),rdis(cm)'
    do i=1,ndata
        write(linebuff,*) dt*dble(i-1),',',acc(i),',',racc(i),',',rvel(i),',',rdis(i)
        call del_spaces(linebuff)
        write(6,'(a)') trim(linebuff)
    end do

stop

contains

    subroutine FUNC(am,ac,ak,ff,x,y,dxdt,dydt)
        !*****************************************************************
        !Function for Runge-Kutta method
        !*****************************************************************
        real(8),intent(in)::am,ac,ak,ff,x,y
        real(8),intent(out)::dxdt,dydt
        dxdt=y
        dydt=ff/am-ac/am*y-ak/am*x
    end subroutine FUNC


    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces

end program f90_acc_rng