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

program f90_EQSP
    !*****************************************************************
    !Acceleration response spetcrum (by Dr. Yorihiko Ohsaki)
    !*****************************************************************
    use defpi
    implicit none
    integer::i,ndata,nnn
    real(8)::dt,h
    real(8)::accmax,velmax,dismax
    real(8),allocatable::acc(:),vel(:),dis(:)
    real(8),allocatable::tk(:),resacc(:),resvel(:),resdis(:)
    character::fnameR*50,fnameW*50
    character::dummy*50,linebuff*300

    h=0.05   !Damping factor
    nnn=200  !Number of devision of period span
    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 i=1,nnn
        tk(i)=10.0D0**(log10(2.0D0*dt)+(log10(10.0D0)-log10(2.0D0*dt))/dble(nnn)*dble(i-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 ERES(ndata,dt,h,acc,nnn,tk,resacc,resvel,resdis)    !calculation of response spectrum

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

    stop

contains

    subroutine ERES(nn,dt,h,acc,nt,tk,resacc,resvel,resdis)
        !*************************************************************
        !Calculation of response spectrum
        !*************************************************************
        !nn      : Number od time history data of acceleration
        !dt      : Time increment                                  (input data)
        !h       : Damping factor                                  (input data)
        !acc()   : Time history data of acceleration               (input data)
        !nt      : Number of period span devision (i=1 to nt)      (input data)
        !tk()    : Calculated period (input data)                  (input data)
        !resacc(): Acceleration response spectrum                  (output data)
        !resvel(): Velocity response spectrum                      (output data)
        !resdis(): Displacement response spectrum                  (output data)

        integer,intent(in)::nn,nt
        real(8),intent(in)::dt,h
        real(8),intent(in)::acc(1:nn)
        real(8),intent(in)::tk(1:nt)
        real(8),intent(out)::resacc(1:nt),resvel(1:nt),resdis(1:nt)

        integer::i,m
        real(8)::accmax,velmax,dismax
        real(8)::w,w2,hw,wd,wdt,e,cwdt,swdt
        real(8)::ss,cc,s1,c1,s2,c2,s3,c3
        real(8)::A11,A12,A21,A22,B11,B12,B21,B22
        real(8)::dxf,xf,ddym,ddyf,x,dx,ddx

        do i=1,nt
            w=2.0D0*pi/tk(i)
            w2=w*w
            hw=h*w
            wd=w*sqrt(1.0D0-h*h)
            wdt=wd*dt
            e=exp(-hw*dt)
            cwdt=cos(wdt)
            swdt=sin(wdt)
            A11=e*(cwdt+hw*swdt/wd)
            A12=e*swdt/wd
            A21=-e*w2*swdt/wd
            A22=e*(cwdt-hw*swdt/wd)
            ss=-hw*swdt-wd*cwdt
            cc=-hw*cwdt+wd*swdt
            s1=(e*ss+wd)/w2
            c1=(e*cc+hw)/w2
            s2=(e*dt*ss+hw*s1+wd*c1)/w2
            c2=(e*dt*cc+hw*c1-wd*s1)/w2
            s3=dt*s1-s2
            c3=dt*c1-c2
            B11=-s2/wdt
            B12=-s3/wdt
            B21=(hw*s2-wd*c2)/wdt
            B22=(hw*s3-wd*c3)/wdt
            accmax=2.0D0*hw*abs(acc(1))*dt
            velmax=abs(acc(1))*dt
            dismax=0.0
            dxf=-acc(1)*dt
            xf=0.0
            do m=2,nn
                ddym=acc(m)
                ddyf=acc(m-1)
                x=A12*dxf+A11*xf+B12*ddym+B11*ddyf
                dx=A22*dxf+A21*xf+B22*ddym+B21*ddyf
                ddx=-2.0D0*hw*dx-w2*x
                if(accmax<=abs(ddx))accmax=abs(ddx)
                if(velmax<=abs(dx))velmax=abs(dx)
                if(dismax<=abs(x))dismax=abs(x)
                dxf=dx
                xf=x
            end do
            resacc(i)=accmax
            resvel(i)=velmax
            resdis(i)=dismax
        end do
    end subroutine ERES


    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_EQSP
