PROGRAM test

  use, intrinsic :: iso_c_binding 
  IMPLICIT none
  
  INCLUDE 'fftw3.f03'

  real(kind=8), parameter :: pi = acos(-1.0) ! also used from f1()

  ! rest of variables only used from "main" program code
  integer, parameter :: npoints = 65536, ncomplex = npoints / 2 + 1
  real(kind=8), parameter :: start = 0, end = 1
  integer :: i
  real(kind=8) :: t,                                            &
       delta = (end - start) / real(npoints - 1, C_DOUBLE)

  type(C_PTR) :: p_real, p_complex, plan
  
  real(C_DOUBLE), pointer :: arr_real(:)
  complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:)

  character(100) :: arg

  IF (command_argument_count() < 1) THEN
     
     DO i = 1, npoints
        
        t = (i - 1) * delta + start
        
        write(*,*) t, " ", f1(t)
        
     END DO
     
     STOP
  END IF

  call get_command_argument(1, arg)

  IF (trim(arg) /= "dft") THEN

     STOP "wrong argument"
  END IF
  
  p_real = fftw_alloc_real(int(npoints, C_SIZE_T))

  p_complex = fftw_alloc_complex(int(ncomplex, C_SIZE_T))

  call c_f_pointer(p_real,    arr_real,    [npoints])
  call c_f_pointer(p_complex, arr_complex, [ncomplex])
  
  plan = fftw_plan_dft_r2c_1d(int(npoints, C_INT), arr_real,     &
                              arr_complex, FFTW_ESTIMATE)
  
  DO i = 1, npoints

     t = (i - 1) * delta + start
     arr_real(i) = real(f1(t), C_DOUBLE)
     
  END DO

  call fftw_execute_dft_r2c(plan, arr_real, arr_complex)

  DO i = 1, ncomplex

     write(*,*) i - 1, " ", abs(arr_complex(i))
     
  END DO

  call fftw_free(p_real)
  call fftw_free(p_complex)
  
CONTAINS

  PURE real(kind=8) FUNCTION f1(t)
    real(kind=8), intent(in) :: t
    
    f1 = sin(2 * pi * t * 200) + 2 * sin(2 * pi * t * 400)

  END FUNCTION f1
  
END PROGRAM test