diff options
Diffstat (limited to 'src/fourier.f90')
-rw-r--r-- | src/fourier.f90 | 116 |
1 files changed, 116 insertions, 0 deletions
diff --git a/src/fourier.f90 b/src/fourier.f90 new file mode 100644 index 0000000..be0a16f --- /dev/null +++ b/src/fourier.f90 @@ -0,0 +1,116 @@ +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 = .true. + call random_seed() + call random_number(randvals) + 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) /= "transformed")) & + STOP "bad second argument (should be 'dft' or 'transformed')" + + 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 + STOP "not implemented" + 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 |