      program dft1
      
      
      integer nmax, n, half
      common n, half
      parameter (nmax = 100)
      complex x(1: nmax),x2(1: nmax),x3(1:nmax),x4(1:nmax),x5(1:nmax)
      real PI
      complex f(1: nmax),I, f2(1:nmax),f3(1:nmax),f4(1:nmax),f5(1:nmax)
      integer sign, outer, inner,q,k,j,r,u,y,h
      PI = 4*ATAN(1.0)
      half = n / 2.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 10 q = 1, n
         x(q) = q
         write(*,*) x(q)
 10   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 50 j = 1, n
         f2(j) = f(j)
         
         x2(j) = 0.0
 50   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 60 r = 1, n
         x2(r) = x2(r)/n
         write(*,*) x2(r)
 60   continue
      
      do 70 u = 1, half
         x3(u) = x(u * 2)
         x4(u) = x(u * 2 - 1)
 70   continue
      x3(1) = x(2)
      x3(2) = x(4)
      x3(3) = x(6)
      x3(4) = x(8)
      x4(1) = x(1)
      x4(2) = x(3)
      x4(3) = x(5)
      x4(4) = x(7)
      
      write(*,*) x3(1),x3(2),x3(3),x3(4),x4(1),x4(2),x4(3),x4(4)
      write(*,*) x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8)

      write(*,*) ' The results from the FFT are '
      call fft(x3,f3,x4,f4)
      do 90 h = 1, half
         f5(2 * h) = f3(h)
         f5(2 * h - 1) = f4(h)
 90   continue
      
      do 100 y = 1, n
         write(*,*) f5(y)
 100  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 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 fft(array1,array2,array3,array4)
      integer n, half,g
      common n, half
      complex array1(half), array2(half), array3(half), array4(half),I
      real PI
      half = n / 2.0
      PI = 4.0 * ATAN(1.0)
      I = (0,1)
      
      call dft(array1,array2)
      call dft(array3,array4)
      
      do 80 g = 1, half
         val = 2.0 * g * PI / n
         array4(g) = array4(g) * EXP(I*val)
 80   continue
      
      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 20 outer = 1, n
         out(outer) = 0.0
         do 30 inner = 1, n
            value = sign*2.0*PI*(outer-1)*(inner-1)/n
            out(outer) = out(outer) + in(inner)*EXP(I*value)
 30      continue
 20   continue
      return
      end
      
