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


program f90_acc_fft
    !*****************************************************************
    !Acceleration Response Spectrum (frequency domain)
    !*****************************************************************
    use defpi
    implicit none
    integer::ndata,i
    real(8),allocatable::acc(:)     !Acceleration
    real(8),allocatable::racc(:)    !Acceleration
    real(8),allocatable::rvel(:)    !Velocity
    real(8),allocatable::rdis(:)    !Displacement
    real(8)::hh                     !Damping factor
    real(8)::Tp                     !Natural period
    real(8)::dt                     !Time increment
    character::fnameR*50            !Input & output file name
    character::dummy*50,linebuff*300
    real(8)::w0

    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))
    allocate(racc(1:ndata))
    allocate(rvel(1:ndata))
    allocate(rdis(1:ndata))
    do i=1,ndata
        read(11,*) acc(i)
    end do
    close(11)

    w0=2.0D0*pi/Tp
    call ARSFFT(ndata,acc,dt,hh,w0,racc,rvel,rdis)

    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)+acc(i),',',rvel(i),',',rdis(i)
        call del_spaces(linebuff)
        write(6,'(a)') trim(linebuff)
    end do

stop

contains

    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)::w_f,w_a,w_h,w_r

        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)*w(i)-w0*w0,-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),disr(i),disi(i),velr(i),veli(i))
            call Z1timesZ2(-w(i)*w(i),0.0D0,disr(i),disi(i),accr(i),acci(i))
        end do
        open(12,file='_acc_fft.txt',status='replace')
        do i=1,nn/2+1 !Memory for reference output
            w_f=dble(i)/dble(nn)/dt
            w_a=sqrt(ar(i)**2+ai(i)**2)
            w_h=sqrt(tr(i)**2+ti(i)**2)
            w_r=sqrt(accr(i)**2+acci(i)**2)
            write(12,'(4(e15.7))') w_f,w_a,w_h,w_r
        end do
        close(12)
        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
        open(12,file='_acc_fft_ri.txt',status='replace')
        do i=1,ndata
            acc(i)=accr(i)
            vel(i)=veli(i)
            dis(i)=disr(i)
            write(12,'(6(e15.7))') accr(i),acci(i),velr(i),veli(i),disr(i),disi(i)
        end do
        close(12)
    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 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_fft
