program main
    use bremModule
    implicit none

! for the gaussian
    real*8, parameter :: pi2=4.d0*atan(1.d0)
    real*8, parameter :: sig=1.d-5
    real*8, parameter :: const=1.d0/(sig*sqrt(2*pi2))
    real*8 tempReal

! for the FFT
    real*8, parameter :: nmax=3.14*10, dn=1.0d-3
    integer, parameter :: nLength=2*nmax/dn
    complex*16 nListIn(0:nlength),nListOut(0:nLength)
    integer j

    OPEN(UNIT=15, FILE='out/in.txt')
    OPEN(UNIT=16, FILE='out/outReal.txt')
    OPEN(UNIT=17, FILE='out/outImag.txt')

    do j=-nLength/2,nLength/2-1
! true delta function (inverse transform)
       if (j*dn.eq.1.5) then
          nListIn(mod(j+nLength,nLength))=dcmplx(1.d0/dn,0.0d0)
       else
          nListIn(mod(j+nLength,nLength))=(0.0d0,0.0d0)
       end if

! narrow gaussian (backward transform)
!      tempReal=const*exp(-(j*dn-1.5d0)**2/(2.d0*sig**2))
!      nListIn(mod(j+nLength,nLength))=dcmplx(tempReal,0.d0)

! testing the forward transform (constant function)
!       nListIn(mod(j+nLength,nLength))=dcmplx( 1.d0,0.d0 )

    end do

    nListOut=FFTB(nListIn,nLength)
!    nListOut=FFTF(nListIn,nLength)

    do j=-nLength/2,nLength/2-1
       write(15,*) j*dn,real( nListIn (mod(j+nLength,nLength)) )
       write(16,*) j*dn,real( nListOut(mod(j+nLength,nLength)) )
       write(17,*) j*dn,imag( nListOut(mod(j+nLength,nLength)) )
    end do

    close(15)
    close(16)
    close(17)
end
