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


program f90_FSP
    !*****************************************************************
    !Fourier Spectrum
    !*****************************************************************
    use defpi
    implicit none
    integer::i,ndata,nn
    real(8)::dt,df,band
    real(8),allocatable::accinp(:)
    real(8),allocatable::xr(:)
    real(8),allocatable::xi(:)
    real(8),allocatable::fs0(:)
    real(8),allocatable::fs1(:)
    real(8),allocatable::frq(:)
    character::fnameR*50,fnameW*50
    character::strcom*200,dummy*50
    character::linebuff*1000

    call getarg(1,fnameR)
    call getarg(2,fnameW)
    call getarg(3,dummy)
    read(dummy,*) band

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) dummy,dt
    read(11,*) dummy,ndata
    nn=2
    do
        nn=nn*2 !to add zeros if necessary
        if(ndata*1<=nn)exit
    end do
    allocate(accinp(1:nn))
    allocate(xr(1:nn))
    allocate(xi(1:nn))
    allocate(fs0(1:nn/2+1))
    allocate(fs1(1:nn/2+1))
    allocate(frq(1:nn/2+1))
    do i=1,nn  !to set number of data to powers of 2
        accinp(i)=0.0D0
        xr(i)=0.0D0
        xi(i)=0.0D0
    end do
    do i=1,ndata
        read(11,*) accinp(i)
        xr(i-1)=accinp(i)
    end do
    close(11)

    call FFT(nn,xr,xi)  !Fourier transform
    do i=1,nn  !to divide return values by number of data
        xr(i)=xr(i)/dble(nn)
        xi(i)=xi(i)/dble(nn)
    end do

    df=1.0D0/dt/dble(nn)
    do i=1,nn/2+1 !Fourier spectrum and frequency
        fs0(i)=dt*dble(nn)*sqrt(xr(i)**2.0D0+xi(i)**2.0D0)
        fs1(i)=fs0(i)
        frq(i)=dble(i)*df
    end do
    call SWIN(nn,fs1,df,band)  !Smoothing of spectrum

    open (12, file=fnameW, status='replace')
    write(12,'(a)') trim(strcom)
    write(12,'(a)') 'i,frequency(Hz),sp0(gal*sec),sp1(gal*sec)'
    do i=1,nn/2+1
        write(linebuff,*) i,',',frq(i),',',fs0(i),',',fs1(i)
        call del_spaces(linebuff)
        write(12,'(a)') trim(linebuff)
    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 SWIN(nn,fs,df,band)
        !**************************************************
        !Smoothing of Fourier spectrum (Parzen window)
        !**************************************************
        !nn    : number of input acceleration data
        !fs()  : Fourier spectrum
        !df    : frequency interval of spectrum (Hz)
        !band  : width of band (Hz)
        !nfold : number of spectrum data

        integer,intent(in)::nn
        real(8),intent(in)::df,band
        real(8),intent(inout)::fs(1:nn/2+1)

        real(8),allocatable::w(:)
        real(8),allocatable::g(:)
        real(8),allocatable::g1(:)
        real(8),allocatable::g2(:)
        real(8)::tt,udf,difs,dif,s
        integer::i,j,nfold
        integer::lmax,ll,ln,lt,le

        nfold=nn/2+1

        if(band>0.0)then
            tt=1.0D0/df
            udf=280.0D0/151.0D0/band*df
            lmax=int(2.0D0/udf)+1
            if(udf<=0.5D0)then
                ll=lmax*2-1
                ln=ll-1+nfold
                lt=(ll-1)*2+nfold
                le=lt-lmax+1
                allocate(w(1:lmax))
                allocate(g(1:nfold))
                allocate(g1(1:lt))
                allocate(g2(1:le))
                w(1)=0.75D0*udf
                do i=2,lmax
                    dif=0.5D0*pi*dble(i)*udf
                    w(i)=w(1)*(sin(dif)/dif)**4.0D0
                end do
                g(1)=fs(1)*fs(1)/tt
                do i=2,nfold-1
                    g(i)=2.0D0*fs(i)*fs(i)/tt
                end do
                g(nfold)=fs(nfold)*fs(nfold)/tt
                do i=1,lt
                    g1(i)=0.0D0
                end do
                do i=1,nfold
                    g1(ll-1+i)=g(i)
                end do
                do i=lmax,le
                    s=w(1)*g1(i)
                    do j=2,lmax
                        s=s+w(j)*(g1(i-j+1)+g1(i+j-1))
                    end do
                    g2(i)=s
                end do
                do j=2,lmax
                    g2(ll+j-1)=g2(ll+j-1)+g2(ll-j+1)
                    g2(ln-j+1)=g2(ln-j+1)+g2(ln+j-1)
                end do
                do i=1,nfold
                    g(i)=g2(ll-1+i)
                end do
                fs(1)=sqrt(g(1)*tt)
                do i=2,nfold-1
                    fs(i)=sqrt(g(i)*tt/2.0D0)
                end do
                fs(nfold)=sqrt(g(nfold)*tt)
            end if
        end if
    end subroutine SWIN


    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_FSP