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

program f90_DECO
    use defpi
    implicit none
    integer::i,ndata,nn
    real(8)::dt
    real(8),allocatable::accinp1(:)
    real(8),allocatable::xr1(:)
    real(8),allocatable::xi1(:)
    real(8),allocatable::accinp2(:)
    real(8),allocatable::xr2(:)
    real(8),allocatable::xi2(:)
    real(8),allocatable::accout(:)
    real(8),allocatable::cr(:)
    real(8),allocatable::ci(:)
    real(8)::tempr,tempi
    character::fnameR1*50,fnameR2*50,fnameW*50
    character::strcom*200,dummy*50
    character::linebuff*1000

    call getarg(1,fnameR1)
    call getarg(2,fnameR2)
    call getarg(3,fnameW)

    open(11,file=fnameR1,status='old')
    read(11,'(a)') strcom
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    nn=2
    do
        nn=nn*2 !必要に応じて後続の0を追加
        if(ndata*1<=nn)exit
    end do
    allocate(accinp1(0:nn-1))
    allocate(accinp2(0:nn-1))
    allocate(accout(0:nn-1))
    allocate(xr1(0:nn-1))
    allocate(xi1(0:nn-1))
    allocate(xr2(0:nn-1))
    allocate(xi2(0:nn-1))
    allocate(cr(0:nn-1))
    allocate(ci(0:nn-1))
    do i=1,nn  !データ数を2の累乗にセット
        accinp1(i-1)=0.0D0
        xr1(i-1)=0.0D0
        xi1(i-1)=0.0D0
        accinp2(i-1)=0.0D0
        xr2(i-1)=0.0D0
        xi2(i-1)=0.0D0
    end do
    do i=1,ndata
        read(11,*) accinp1(i-1)
        xr1(i-1)=accinp1(i-1)
    end do
    close(11)

    open(11,file=fnameR2,status='old')
    read(11,'(a)') strcom
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    do i=1,ndata
        read(11,*) accinp2(i-1)
        xr2(i-1)=accinp2(i-1)
    end do
    close(11)

    call FFT(nn,xr1,xi1)  !Fourier変換
    do i=1,nn  !戻り値をデータ数で除す
        xr1(i-1)=xr1(i-1)/dble(nn)
        xi1(i-1)=xi1(i-1)/dble(nn)
    end do
    call FFT(nn,xr2,xi2)  !Fourier変換
    do i=1,nn  !戻り値をデータ数で除す
        xr2(i-1)=xr2(i-1)/dble(nn)
        xi2(i-1)=xi2(i-1)/dble(nn)
    end do

    do i=1,nn !複素数演算
        call Z1timesZ2(xr1(i-1),xi1(i-1),xr1(i-1),xi1(i-1),tempr,tempi)
        call Z1overZ2(tempr,tempi,xr2(i-1),xi2(i-1),cr(i-1),ci(i-1))
    end do
    do i=1,nn  !フーリエ逆変換用に虚数部の符号を反転
        cr(i-1)=cr(i-1)
        ci(i-1)=-ci(i-1)
    end do
    Call FFT(nn,cr,ci)  !Fourier逆変換

    open (12, file=fnameW, status='replace')
    write(12,'(a)') 'result_acc'
    write(linebuff,*) 'dt,',dt
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    write(linebuff,*) 'ndata,',ndata
    call del_spaces(linebuff)
    write(12,'(a)') trim(linebuff)
    do i=1,ndata
        write(linebuff,*) cr(i-1)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    end do
    close(12)

stop

contains

    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(0:nn-1)
        real(8),intent(inout)::xi(0:nn-1)

        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(0:n/2))
        allocate(c(0:n/2))
        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)=sin(a)
            c(i)=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)-xr(j)
                    b=xi(i)-xi(j)
                    xr(i)=xr(i)+xr(j)
                    xi(i)=xi(i)+xi(j)
                    if(p==0)then
                        xr(j)=a
                        xi(j)=b
                    else
                        xr(j)=a*c(p)+b*s(p)
                        xi(j)=b*c(p)-a*s(p)
                    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)
                xr(i)=xr(j)
                xr(j)=xd
                xd=xi(i)
                xi(i)=xi(j)
                xi(j)=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_DECO