! compile with 
! f95 -r8 -f77 -O4 -o x1 x1.f90 -L/usr/local/lib -lfftw -lrfftw
!
module constants

  complex, parameter :: zi = ( 0.0 , 1.0 )
  real, parameter :: pi = 3.1415926535897931
  
end module constants

module fft

  implicit none

! 
! For scientific applications that require a fast Fourier transform,
! one wishes to use an optimized set of subroutines to obtain 
! good performance.  The best portable and well-documented FFT 
! library that I am familiar with is called "FFTW", and is the 
! one that we will use below.  FFTW is freely available at 
! http://www.fftw.org, where there is also extensive documentation. 
! 
! There is one important usage note, which I include for future
! reference: : If you wish to perform one-dimensional 
! real-to-complex or complex-to-real transforms, I recommend
! using the rfftwnd set of calls, with the number of dimensions set
! to one.  Otherwise, you will likely find that FFTW's ordering of 
! array elements is inconvenient and confusing.  
!
! The FFTW routines are not written in Fortran.  To call them 
! from Fortran programs requires a careful look at the documentation, 
! or you can follow the examples below.  The interface is reasonably
! intuitive.
!

  integer fftw_forward,fftw_backward
  parameter (fftw_forward=-1,fftw_backward=1)
  
  integer fftw_real_to_complex,fftw_complex_to_real
  parameter (fftw_real_to_complex=-1,fftw_complex_to_real=1)
  
  integer fftw_estimate,fftw_measure
  parameter (fftw_estimate=0,fftw_measure=1)
  
  integer fftw_out_of_place,fftw_in_place
  parameter (fftw_out_of_place=0)
  parameter (fftw_in_place=8)

end module fft

module output

contains
  
  subroutine x_out (x, u)
    real, dimension(:), intent (in) :: x, u
    integer :: n, i

    n=size(u)

    open (unit=22,file='x1.x') 
    write(22,*) '#  x                      u'
    do i=1,n
       write(22,*) x(i), u(i)
    end do

    ! use periodicity to make figures pretty
    write(22,*) x(n)+x(2)-x(1), u(1)

    close (22)
  end subroutine x_out

  subroutine k_out (k, a)
    complex, dimension(:), intent (in) :: a
    real, dimension(:), intent (in) :: k
    integer :: n, i

    n = size(a)

    open(unit=21,file='x1.k')
    write(21,*) '#  kx                    a_k'
    do i=1,n
       write(21,*) k(i), real(a(i)), aimag(a(i))
    end do
    close (21)

  end subroutine k_out


end module output

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!   Main Program starts here   !!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program x1

!
! This is how one gets access to variables declared in a module
!
  use fft  
  use output
  use constants

  implicit none

!
! Declare the variables
!

  real :: dx, c, s
  integer :: n, i, kc, ks, j

  complex, dimension(:), allocatable :: out
  real, dimension(:), allocatable :: in, backout, x, kx

  integer plan
  namelist /input/ n, c, s, kc, ks

! 
! Get the input values. 
! 
  open (unit=16,file='x1.in',status='old')
  read (16, input)
  close (16)

!
! Allocate memory.
!
  allocate (in(n), backout(n))
  allocate (out(n/2+1))
  allocate (x(n), kx(n))

!
! Set up x and kx arrays for later plotting
!
  dx = 2.*pi/n
  x(1) = 0.
  do i = 2, n
     x(i) = x(i-1) + dx
  end do

  do i = 1, n/2+1
     kx(i) = real(i-1)
  end do
  
!
! This is the function we will be plotting.  Note the use of 
! the array syntax. 
!
  do j=1,n
     in(j)=c*cos(kc*x(j))+s*sin(ks*x(j))
  end do
!
! initialize FFTW 
!
  call rfftwnd_f77_create_plan(plan,1,n,fftw_forward,fftw_estimate)
!
! Perform FFT
!
  call rfftwnd_f77_one_real_to_complex(plan,in,out)
!
! Delete the contents of 'plan' because we are going to redefine
! it shortly.
!
  call rfftwnd_f77_destroy_plan(plan)

! Normalization by number of elements is left to the user
  out = out/n

  call x_out (x, in)
  call k_out (kx, out)

! It is easy to do the inverse fft because the parameters are set.
! just replace fftw_forward with fftw_backward in the plan call.

  call rfftwnd_f77_create_plan(plan,1,n,fftw_backward,fftw_estimate)
  call rfftwnd_f77_one_complex_to_real(plan,out,backout)
  call rfftwnd_f77_destroy_plan(plan)

  call x_out(x, backout)

end program x1