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


program f90_wave_phase
    !*********************************************
    ! Fourier Phase Spectrum
    !*********************************************
    use defpi
    implicit none
    integer::i,k                       !
    integer::ndata                   !Number of original data
    integer::nn                      !Number of calculated data (powers of 2)
    real(8)::dt                      !Time increment
    real(8),allocatable::accinp(:)   !Inputted original data
    real(8),allocatable::xr(:)       !Real part of calculated data
    real(8),allocatable::xi(:)       !Imaginary part of calculated data
    real(8),allocatable::fp(:)       !Fourier phase spectrum
    real(8),allocatable::frq(:)       !Fourier phase spectrum
    real(8)::pwave
    character::fnameR*50,fnameW*50,fnameP*50
    character::strcom*200,dummy*50,fmt1*200
    character::linebuff*1000

    call getarg(1,fnameR) !Input data file name
    call getarg(2,fnameW) !Phase angle output file name
    call getarg(3,fnameP) !Phase wave output file name

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    !To set number of data
    ! (nn is equal to 2 to the power of positive integer number 
    !  and nn is greater than or equal to number of data)
    nn=2
    do
        nn=nn*2
        if(ndata*1<=nn)exit
    end do
    allocate(accinp(1:nn))
    allocate(xr(1:nn))
    allocate(xi(1:nn))
    allocate(fp(1:nn))
    allocate(frq(1:nn/2+1))
    do k=1,nn/2+1 !set frequency
        frq(k)=1.0D0/dt/dble(nn)*dble(k-1)
    end do
    do i=1,nn
        accinp(i)=0.0D0
        xr(i)=0.0D0
        xi(i)=0.0D0
    end do
    do i=1,ndata
        read(11,*) accinp(i)
        xr(i)=accinp(i)
    end do
    close(11)

    call FFT(nn,xr,xi)  !Fourier transform
    do i=1,nn  !To divide returned values by number of data
        xr(i)=xr(i)/dble(nn)
        xi(i)=xi(i)/dble(nn)
        fp(i)=0.0D0
        if(1.0D-30<abs(xr(i)))then
            fp(i)=atan(xi(i)/xr(i))/pi*180.0D0
            if(xr(i)<0.0D0.and.xi(i)<0.0D0)fp(i)=fp(i)-180.0D0
            if(xr(i)<0.0D0.and.0.0D0<xi(i))fp(i)=180.0D0+fp(i)
        else
            if(0.0D0<xi(i).and.0.0D0<xr(i))fp(i)=90.0D0
            if(xi(i)<0.0D0.and.xr(i)<0.0D0)fp(i)=-90.0D0
            if(0.0D0<xi(i).and.xr(i)<0.0D0)fp(i)=90.0D0
            if(xi(i)<0.0D0.and.0.0D0<xr(i))fp(i)=-90.0D0
        end if
    end do

    fmt1="(i5,2(',',e15.7))"
    open (12, file=fnameW, status='replace')
    write(12,'(a)') 'k,frq(k),fp(k)'

    do k=1,nn/2+1
        write(linebuff,FMT=fmt1) k,frq(k),fp(k)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    end do
    close(12)

    open (12, file=fnameP, status='replace')
    write(12,'(a)') 'Phase wave'
    write(12,*) 'dt,',dt
    write(12,*) 'ndata,',nn
    do i=1,nn
        pwave=0.0D0
        do k=2,nn/2+1
            pwave=pwave+cos(2.0D0*pi*frq(k)*dble(i)*dt+fp(k)/180.0D0*pi)
        end do
        write(12,*) pwave
    end do
    close(12)

stop

contains

    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_wave_phase
