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


program f90_ARS
    !*****************************************************************
    !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::vel(:)       !Velocity
    real(8),allocatable::dis(:)       !Displacement
    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))
    allocate(vel(1:ndata))
    allocate(dis(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 IACC(ndata,dt,acc,vel,dis,accmax,velmax,dismax)     !Integral of acceleration
    call CRAC(ndata,dt,acc,vel,dis,accmax,velmax,dismax)     !Amendment of base line
    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 NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,dt,ff)
        !*****************************************************************
        !Newmark Beta Method
        !*****************************************************************
        real(8),intent(out)::acc2,vel2,dis2
        real(8),intent(in)::acc1,vel1,dis1
        real(8),intent(in)::am,ac,ak,dt,ff
!        real(8)::beta=1.0D0/4.0D0  !Average acceleration method
        real(8)::beta=1.0D0/6.0D0  !Linear acceleration method

        acc2=(ff-ac*(vel1+0.5D0*dt*acc1)-ak*(dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1))&
             /(am+ac*0.5D0*dt+ak*beta*dt*dt)
        vel2=vel1+0.5D0*dt*(acc1+acc2)
        dis2=dis1+dt*vel1+(0.5D0-beta)*dt*dt*acc1+beta*dt*dt*acc2
    end subroutine NEWMARK


    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,md
        real(8)::acc2,vel2,dis2
        real(8)::acc1,vel1,dis1
        real(8)::am,ac,ak,f1,f2,fm,ddt
        real(8)::racc,rvel,rdis
        real(8)::acmax,vemax,dimax

        md=5
        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

            acc1=0.0D0
            vel1=0.0D0
            dis1=0.0D0
            acmax=acc1
            vemax=vel1
            dimax=dis1
            do i=1,ndata
                f1=0.0D0
                if(1<=i)f1=-am*acc(i-1)
                f2=-am*acc(i)
                do j=1,md
                    fm=f1+(f2-f1)/dble(md)*dble(j)
                    call NEWMARK(acc2,vel2,dis2,acc1,vel1,dis1,am,ac,ak,ddt,fm)
                    acc1=acc2
                    vel1=vel2
                    dis1=dis2
                end do
                racc=acc2+acc(i)     !Absolute acceleration
                rvel=vel2            !Relative velocity
                rdis=dis2            !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 IACC(nn,dt,acc,vel,dis,accmax,velmax,dismax)
        !*****************************************************************
        !Integral of acceleration
        !*****************************************************************
        !nn    : Number of time history data of acceleration
        !dt    : Time increment                    (input data)
        !acc() : Time history data of Acceleration (input data)
        !vel() : Time history data of Velocity     (output data)
        !dis() : Time history data of Displacement (output data)
        !accmax: Maximum value of acceleration     (output data)
        !velmax: Maximum value of velocity         (output data)
        !dismax: Maximum value of displacement     (output data)

        integer,intent(in)::nn
        real(8),intent(in)::dt
        real(8),intent(in)::acc(1:nn)
        real(8),intent(out)::vel(1:nn),dis(1:nn)
        real(8),intent(out)::accmax,velmax,dismax

        integer::i

        accmax=0.0
        velmax=0.0
        dismax=0.0
        vel(1)=0.0
        dis(1)=0.0
        do i=2,nn
            vel(i)=vel(i-1)+(acc(i-1)+acc(i))*dt/2.0D0
            dis(i)=dis(i-1)+vel(i-1)*dt+(acc(i-1)/3.0D0+acc(i)/6.0D0)*dt*dt
            if(accmax<abs(acc(i)))accmax=abs(acc(i))
            if(velmax<abs(vel(i)))velmax=abs(vel(i))
            if(dismax<abs(dis(i)))dismax=abs(dis(i))
        end do
    end subroutine IACC


    subroutine CRAC(nn,dt,acc,vel,dis,accmax,velmax,dismax)
        !*****************************************************************
        !Amendment of base line
        !*****************************************************************
        !nn    : Number of time history data of acceleration
        !dt    : Time increment                    (input data)
        !acc() : Time histori data of acceleration (input & output data)
        !vel() : Time history of velocity          (input data)
        !dis() : Time history data of displacement (input & output data)
        !accmax: Maximum value of acceleration     (input data)
        !velmax: Maximum value of velocity         (input data)
        !dismax: Maximum value of displacement     (input data)

        integer,intent(in)::nn
        real(8),intent(in)::dt,accmax,velmax,dismax
        real(8),intent(in)::vel(1:nn)
        real(8),intent(inout)::dis(1:nn)
        real(8),intent(inout)::acc(1:nn)

        integer::i
        real(8)::tt,t,sum,a0,a1,acmax,coef

        tt=dble(nn-1)*dt
        t=0.0D0
        do i=1,nn
            dis(i)=dis(i)*(3.0D0*tt-2.0D0*t)*t*t
            t=t+dt
        end do
        sum=(dis(1)+dis(nn))/2.0D0
        do i=2,nn-1
            sum=sum+dis(i)
        end do
        sum=sum*dt
        a1=28.0D0/13.0D0/tt/tt*(2.0D0*vel(nn)-15.0D0/tt**5.0D0*sum)
        a0=vel(nn)/tt-a1/2.0D0*tt
        t=0.0D0
        acmax=0.0D0
        do i=1,nn
            acc(i)=acc(i)-a0-a1*t
            if(acmax<abs(acc(i)))acmax=abs(acc(i))
            t=t+dt
        end do
        coef=accmax/acmax
        do i=1,nn
            acc(i)=acc(i)*coef
        end do
    end subroutine CRAC


    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