aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/f1.f9081
-rw-r--r--src/fourier.f90116
2 files changed, 116 insertions, 81 deletions
diff --git a/src/f1.f90 b/src/f1.f90
deleted file mode 100644
index 7f07afa..0000000
--- a/src/f1.f90
+++ /dev/null
@@ -1,81 +0,0 @@
-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
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