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


program f90_ars_rng
    !*****************************************************************
    !Acceleration Response Spectrum (time domein)
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i,k,nnn
    real(8)::accmax,velmax,dismax     !Maximum values of acceleration, velocity and displacement
    real(8),allocatable::acc(:)       !Acceleration (input data)
    real(8),allocatable::resacc(:)    !Acceleration
    real(8),allocatable::resvel(:)    !Velocity
    real(8),allocatable::resdis(:)    !Displacement
    real(8),allocatable::Tk(:)        !Natural period
    real(8)::hh                       !Damping factor
    real(8)::dt                       !Time increment
    character::fnameR*50,fnameW*50    !Input & output file name
    character::dummy*50,linebuff*300

    nnn=200
    hh=0.05

    allocate(Tk(1:nnn))
    allocate(resacc(1:nnn))
    allocate(resvel(1:nnn))
    allocate(resdis(1:nnn))

    call getarg(1,fnameR)  !Input data file name (time history data of acceleration)
    call getarg(2,fnameW)  !Output file name

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

    do k=1,nnn
        Tk(k)=10.0D0**(log10(2.0*dt)+(log10(10.0)-log10(2.0*dt))/dble(nnn)*dble(k-1))
    end do
    call CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis) !Calculation of response

    open (12, file=fnameW, status='replace')
    write(12,'(a)') trim(linebuff)
    write(linebuff,*) 'nt,',nnn
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    write(12,'(a)') 'period(sec),acc(gal),vel(kine),dis(cm)'
    do k=1,nnn
        write(linebuff,*) Tk(k),',',resacc(k),',',resvel(k),',',resdis(k)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    end do
    close(12)

    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 CALRES(ndata,nnn,dt,hh,acc,Tk,resacc,resvel,resdis)
        !*****************************************************************
        !Calculation of response
        !*****************************************************************
        integer,intent(in)::ndata,nnn
        real(8),intent(in)::dt,hh
        real(8),intent(in)::acc(1:ndata),Tk(1:nnn)
        real(8),intent(out)::resacc(1:nnn),resvel(1:nnn),resdis(1:nnn)

        integer::i,j,k
        real(8)::am,ac,ak,f1,f2,fm1,fm2,fmm
        real(8)::racc,rvel,rdis
        real(8)::acmax,vemax,dimax
        real(8)::x1,y1,x2,y2,ddt
        real(8)::kx1,kx2,kx3,kx4
        real(8)::ky1,ky2,ky3,ky4

        ddt=dt/dble(md)
        do k=1,nnn
            am=1.0D0                       !Mass
            ak=4.0D0*pi*pi*am/Tk(k)/Tk(k)  !Stiffness
            ac=2.0D0*hh*sqrt(ak*am)        !Damping constant

            racc=0.0D0
            rvel=0.0D0
            rdis=0.0D0
            acmax=racc
            vemax=rvel
            dimax=rdis
            do i=1,ndata
                f1=0.0D0
                if(2<=i)f1=-am*acc(i-1)
                f2=-am*acc(i)
                x1=rdis
                y1=rvel
                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
                racc=f2/am-ac/am*y2-ak/am*x2
                racc=racc+acc(i)   !Absolute acceleration
                rvel=y2            !Relative velocity
                rdis=x2            !Relative displacement
                if(acmax<abs(racc))acmax=abs(racc)
                if(vemax<abs(rvel))vemax=abs(rvel)
                if(dimax<abs(rdis))dimax=abs(rdis)
            end do
            resacc(k)=acmax
            resvel(k)=vemax
            resdis(k)=dimax
        end do
    end subroutine CALRES


    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_ars_rng

