      program dft1

      
      integer nmax, n
      common n
      parameter (nmax = 100)
      complex x(1: nmax),x2(1: nmax)
      real PI
      complex f(1: nmax),I, f2(1:nmax)
      integer sign, outer, inner,j,k,l,r
      PI = 4*ATAN(1.0)
      
      write(*,*) 'What is the number of data points?'
      read(*,*) n
      if (n .gt. nmax) stop 'That size is not allowed. '
      write(*,*) 'The input values to be calculated on the '
      write(*,*) 'Fourier transform are: x(j) = ' , 'Pi = ', PI
      do 30 j = 1, n
         x(j) = j
         write(*,*) x(j)
 30   continue
      
      call dft(x,f)
      
      write(*,*) 'The results from the Fourier transform are f(k) = '
      do 40 k = 1, n
         write(*,*) f(k)
 40   continue


      do 60 j = 1, n
         f2(j) = f(j)

         x2(j) = 0.0
 60   continue
      
      write(*,*) 'The input values for the Inverse Fourier transform '
      write(*,*) 'are the output values f from the last calculation'
      call dftinverse(f2,x2)
      
      write(*,*) 'The results from the Inverse Fourier transform are '
      write(*,*) ' x2(r) = '
      do 50 r = 1, n
         x2(r) = x2(r)/n
         write(*,*) x2(r)
 50   continue


      return
      end

**********************************************************************
*     subroutine dft(input, output) this subroutine calculates the 
*     generalized discrete fourier transform, using dftbody,where sign
*     is equal to 1 
*     usage: dft(input, output)
*     output: out(n), gives an array of size n
**********************************************************************
      
      subroutine dft(array1,array2)
      integer n
      common n
      call dftbody(array1,array2,1)
      return 
      end


**********************************************************************
*     subroutine dftinverse(input, output) this subroutine calculates the 
*     generalized discrete fourier transform, using dftbody,where sign
*     is equal to 1 
*     usage: dft(input, output)
*     output: out(n), gives an array of size n
**********************************************************************
      
      subroutine dftinverse(array1 , array2)
      integer n,d
      common n
      call dftbody(array1,array2,-1)

      return 
      end


**********************************************************************
*     subroutine dftbody(in, new, sign) this subroutine calculates the 
*     generalized discrete fourier transform where sign is either 1 or
*     -1 depending on whether you are calculating the fourier or the 
*     inverse fourier transform, respectively
*     usage: call dftbody(in, out2, sign)
*     output: out(n), gives an array of size n
**********************************************************************
      
      subroutine dftbody(in,out,sign)
      integer n
      common n
      real PI
      complex in(1:n)
      complex out(1:n),I
      integer outer, inner,sign
      PI = 4.0*ATAN(1.0)
      I = (0,1)
      do 10 outer = 1, n
         out(outer) = 0.0
         do 20 inner = 1, n
            value = sign*2.0*PI*(outer-1)*(inner-1)/n
            out(outer) = out(outer) + in(inner)*EXP(I*value)
 20      continue
 10   continue
      return
      end


