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)

  real(kind=8) :: val, randvals(npoints)
  
  type(C_PTR) :: p_real, p_complex, plan_f, plan_b ! forward,backward
  
  real(C_DOUBLE), pointer :: arr_real(:)
  complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:)

  character(100) :: arg

  logical :: which_fun
  
  IF (command_argument_count() < 1) STOP "argument required"

  call get_command_argument(1, arg)

  IF      (arg == "f1") THEN
     which_fun = .true.
  ELSE IF (arg == "f2") THEN
     which_fun = .false.
     call random_seed()
     call random_number(randvals)
     randvals = randvals / 50
  ELSE
     STOP "bad first argument (should be 'f1' or 'f2')"
  END IF
  
  IF (command_argument_count() < 2) THEN
     DO i = 1, npoints
        
        t = (i - 1) * delta + start
        
        IF (which_fun) val = f1(t)   ! f1
        
        IF (.not. which_fun) THEN    ! f2
           val = f2(t) + randvals(i)
        END IF
        
        write(*,*) t, " ", val
        
     END DO
     STOP
  END IF

  call get_command_argument(2, arg)
  
  IF ((trim(arg) /= "dft") .and. (trim(arg) /= "filtered"))      &
       STOP "bad second argument (should be 'dft' or 'filtered')"

  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_f = fftw_plan_dft_r2c_1d(int(npoints, C_INT), arr_real,   &
                                arr_complex, FFTW_ESTIMATE)     
  plan_b = fftw_plan_dft_c2r_1d(int(npoints, C_INT), arr_complex,&
                                arr_real, FFTW_ESTIMATE)
          
  IF (which_fun) THEN     ! f1
     arr_real = real(f1((/(i, i = 1, npoints, 1)/) * delta + start),&
                     C_DOUBLE)
  ELSE                    ! f2
     arr_real = real(f2((/(i, i = 1, npoints, 1)/) * delta + start  &
                        + randvals), C_DOUBLE)
  END IF
     
  call fftw_execute_dft_r2c(plan_f, arr_real, arr_complex)

  IF (trim(arg) == "dft") THEN
     DO i = 1, ncomplex
        
        write(*,*) i - 1, " ", abs(arr_complex(i))
     END DO
  ELSE
     DO i = 1, ncomplex
        
        IF (abs(arr_complex(i)) < 50) arr_complex(i) = 0
     END DO

     call fftw_execute_dft_c2r(plan_b, arr_complex, arr_real)

     arr_real = arr_real / npoints      ! normalize
     
     DO i = 1, npoints
        
        write(*,*) (i - 1) * delta + start, " ", arr_real(i)
     END DO

  END IF
     
  call fftw_free(p_real)
  call fftw_free(p_complex)

CONTAINS

  PURE ELEMENTAL 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

  
  PURE ELEMENTAL real(kind=8) FUNCTION f2(t)
    real(kind=8), intent(in) :: t
    
    f2 = cos(2 * pi * t)

  END FUNCTION f2
  
END PROGRAM test