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


program f90_ARSFFT
    !*****************************************************************
    !Acceleration Response Spectrum (frequency domain)
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i,k,nnn
    real(8)::accmax,velmax,dismax   !Maximum values
    real(8),allocatable::acc(:)     !Acceleration
    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 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 CALRES(ndata,nnn,dt,hh,accinp,Tk,resacc,resvel,resdis)
        integer,intent(in)::ndata,nnn
        real(8),intent(in)::dt,hh
        real(8),intent(in)::accinp(1:ndata),Tk(1:nnn)
        real(8),intent(out)::resacc(1:nnn),resvel(1:nnn),resdis(1:nnn)

        integer::i,k
        real(8)::w0,acmax,vemax,dimax
        real(8)::acc(1:ndata),vel(1:ndata),dis(1:ndata)

        do k=1,nnn
            w0=2.0D0*pi/Tk(k)
            call ARSFFT(ndata,accinp,dt,hh,w0,acc,vel,dis)
            acmax=0.0D0
            vemax=0.0D0
            dimax=0.0D0
            do i=1,ndata
                if(acmax<abs(acc(i)+accinp(i)))acmax=abs(acc(i)+accinp(i))
                if(vemax<abs(vel(i)))vemax=abs(vel(i))
                if(dimax<abs(dis(i)))dimax=abs(dis(i))
            end do
            resacc(k)=acmax
            resvel(k)=vemax
            resdis(k)=dimax
        end do
    end subroutine CALRES


    subroutine ARSFFT(ndata,accinp,dt,hh,w0,acc,vel,dis)
        integer,intent(in)::ndata
        real(8),intent(in)::accinp(1:ndata)
        real(8),intent(in)::dt,hh,w0
        real(8),intent(out)::acc(1:ndata)
        real(8),intent(out)::vel(1:ndata)
        real(8),intent(out)::dis(1:ndata)

        integer::nn,i
        real(8),allocatable::ar(:),ai(:)
        real(8),allocatable::tr(:),ti(:)
        real(8),allocatable::accr(:),acci(:)
        real(8),allocatable::velr(:),veli(:)
        real(8),allocatable::disr(:),disi(:)
        real(8),allocatable::w(:)
        real(8)::fk
        real(8)::tempr,tempi

        nn=2
        do
            nn=nn*2 !set data number to a power of two (data number must be a power oftwo)
            if(ndata*1<=nn)exit
        end do
        allocate(ar(1:nn))
        allocate(ai(1:nn))
        allocate(tr(1:nn))
        allocate(ti(1:nn))
        allocate(accr(1:nn))
        allocate(acci(1:nn))
        allocate(velr(1:nn))
        allocate(veli(1:nn))
        allocate(disr(1:nn))
        allocate(disi(1:nn))
        allocate(w(1:nn))
        do i=1,nn
            ar(i)=0.0D0
            ai(i)=0.0D0
            tr(i)=0.0D0
            ti(i)=0.0D0
        end do
        do i=1,ndata
            ar(i)=accinp(i)
        end do
        call FFT(nn,ar,ai)  !Fourier transform
        do i=1,nn  !To divide returned value by data number
            ar(i)=ar(i)/dble(nn)
            ai(i)=ai(i)/dble(nn)
        end do
        do i=1,nn/2+1 !frequency response function
            fk=dble(i)/dble(nn)/dt
            w(i)=2.0D0*pi*fk
            call Z1overZ2(1.0D0,0.0D0,w(i)**2.0D0-w0**2.0D0,-2.0D0*hh*w0*w(i),tr(i),ti(i))
        end do
        do i=2,nn/2+1
            w(nn-i)=w(i)
            tr(nn-i)= tr(i)
            ti(nn-i)=-ti(i)
        end do
        do i=1,nn !Fourier transform of displacement, velocity and acceleration
            call Z1timesZ2(ar(i),ai(i),tr(i),ti(i),disr(i),disi(i))
            call Z1timesZ2(0.0D0,w(i),ar(i),ai(i),tempr,tempi)
            call Z1timesZ2(tempr,tempi,tr(i),ti(i),velr(i),veli(i))
            call Z1timesZ2(-w(i)*w(i),0.0D0,ar(i),ai(i),tempr,tempi)
            call Z1timesZ2(tempr,tempi,tr(i),ti(i),accr(i),acci(i))
        end do
        do i=1,nn  !Inversion of sign of imaginary parts
            disi(i)=-disi(i)
            veli(i)=-veli(i)
            acci(i)=-acci(i)
        end do
        Call FFT(nn,disr,disi)  !Inverse Fourier transform
        Call FFT(nn,velr,veli)  !Inverse Fourier transform
        Call FFT(nn,accr,acci)  !Inverse Fourier transform
        do i=1,ndata
            acc(i)=accr(i)
            vel(i)=velr(i)
            dis(i)=disr(i)
        end do
    end subroutine ARSFFT


    subroutine Z1timesZ2(z1r,z1i,z2r,z2i,zzr,zzi)
        real(8),intent(in)::z1r,z1i,z2r,z2i
        real(8),intent(out)::zzr,zzi

        zzr=z1r*z2r-z1i*z2i
        zzi=z1r*z2i+z1i*z2r
    end subroutine Z1timesZ2


    subroutine Z1overZ2(z1r,z1i,z2r,z2i,zzr,zzi)
        real(8),intent(in)::z1r,z1i,z2r,z2i
        real(8),intent(out)::zzr,zzi

        zzr=(z1r*z2r+z1i*z2i)/(z2r*z2r+z2i*z2i)
        zzi=(z1i*z2r-z1r*z2i)/(z2r*z2r+z2i*z2i)
    end subroutine Z1overZ2



    subroutine FFT(nn,xr,xi)
    !**************************************************
    !Fast Fourier Transform, Inverse transform
    !**************************************************
    !nn   : Number of data (powers of 2)
    !xr() : Real part of i/o data
    !xi() : Imaginary part of i/o data
        integer,intent(in)::nn
        real(8),intent(inout)::xr(1:nn)
        real(8),intent(inout)::xi(1:nn)

        integer::g,h,i,j,k
        integer::l,m,n,p,q
        real(8)::a,b,xd
        real(8),allocatable::s(:)
        real(8),allocatable::c(:)

        n=nn
        allocate(s(1:n/2+1))
        allocate(c(1:n/2+1))
        i=0;j=0;k=0;l=0;p=0;h=0;g=0;q=0
        m=int(log(dble(n))/log(2.0D0)+1.0D0)
        a=0.0D0
        b=pi*2.0D0/dble(n)
        do i=0,n/2
            s(i+1)=sin(a)
            c(i+1)=cos(a)
            a=a+b
        end do
        l=n
        h=1
        do g=1,m
            l=l/2
            k=0
            do q=1,h
                p=0
                do i=k,l+k-1
                    j=i+l
                    a=xr(i+1)-xr(j+1)
                    b=xi(i+1)-xi(j+1)
                    xr(i+1)=xr(i+1)+xr(j+1)
                    xi(i+1)=xi(i+1)+xi(j+1)
                    if(p==0)then
                        xr(j+1)=a
                        xi(j+1)=b
                    else
                        xr(j+1)=a*c(p+1)+b*s(p+1)
                        xi(j+1)=b*c(p+1)-a*s(p+1)
                    end if
                    p=p+h
                end do
                k=k+l+l
            end do
            h=h+h
        end do
        j=n/2
        do i=1,n-1
            k=n
            if(j<i)then
                xd=xr(i+1)
                xr(i+1)=xr(j+1)
                xr(j+1)=xd
                xd=xi(i+1)
                xi(i+1)=xi(j+1)
                xi(j+1)=xd
            end if
            k=k/2
            do while(j>=k)
                j=j-k
                k=k/2
                if(k==0)exit
            end do
            j=j+k
        end do
    end subroutine FFT


    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_ARSFFT