PROGRAM fourier

  implicit none

  INTEGER, PARAMETER ::  n  = 20
  INTEGER, PARAMETER ::  a1 = 8
  INTEGER, PARAMETER ::  a2 = 5
  INTEGER, PARAMETER ::  b1 = 4
  INTEGER, PARAMETER ::  b2 = 2

  double precision :: pi= acos(-1.d0)

  complex(kind=4),  DIMENSION(n+1) :: usmall
  COMPLEX(kind=4),  dimension(n+1) :: ubig, uinverse

  integer, parameter :: isign01 = 1
  integer, parameter :: isign02 = -1
  integer     :: k


 WRITE(6,*) 'This is fourier_Program'

 !discrete function (1)
 DO k=0,n-1
   usmall(k)=cmplx(a1*cos((pi*2*b1*(k))/n)+a2*sin(2*pi*b2*(k)/n), 0)
 end do

    !aufgabe a/c
 call slowfourier(usmall, isign01, ubig)

  write(*,*)
  write(*,*)
  write(*,*) "Now inverse Transform"
  write(*,*)

 !exercise d
 call slowfourier(ubig/n, isign02, uinverse)


 WRITE(6,*) 'fourier_Program: DONE.'


END PROGRAM



subroutine slowfourier(usub, isign, blubb)
  IMPLICIT NONE

  INTEGER, PARAMETER ::  n1 = 20
  double precision :: pi= acos(-1.d0)

  INTEGER :: ncid

  INTEGER :: nstep1, k1
  integer, intent(in):: isign
  complex(kind=4),  DIMENSION(n1+1), intent(in) :: usub
  complex(kind=4) :: ufct
  complex(kind=4), dimension(n1+1) :: blubb
  !complex(kind=4) :: ni = cmplx(n,n)





    do nstep1=1,n1

    ufct = cmplx(0,0)

    do k1=0, n1-1
      ufct=ufct+ usub(k1)*cmplx(cos(isign*2*pi*k1*(nstep1)/n1), sin(isign*2*pi*k1*(nstep1)/n1))
    end do

    !need this for inverse
    blubb(nstep1)= ufct

   ! write output
    write(*,*) "n=", nstep1
    write(*,*) "ufct=(Re,Im)= ", ufct
    write(*,FMT="(A15,2(F8.5))") "ufct/N= ", ufct/n1

   end do


end subroutine