perform fft, print abs() of resulting values

1 files changed, 24 insertions, 9 deletions

diff --git a/src/f1.f90 b/src/f1.f90
index f03b9bf..82e60b1 100644
--- a/src/f1.f90
+++ b/src/f1.f90
@@ -5,31 +5,46 @@ PROGRAM test
   
   INCLUDE 'fftw3.f03'
 
-  integer, parameter :: npoints = 1024
+  integer, parameter :: npoints = 1024, ncomplex = npoints / 2 + 1
   real(C_DOUBLE), parameter :: start = 0, end = 1
   integer :: i
   real(C_DOUBLE) :: t,                                            &
        delta = (end - start) / real(npoints - 1, C_DOUBLE)
 
-  type(C_PTR) :: p_in!, p_out, plan
+  type(C_PTR) :: p_real, p_complex, plan
 
-  real(C_DOUBLE), pointer :: arr_in(:)
-!  complex(C_DOUBLE_COMPLEX), pointer :: out_arr(npoints)
+  real(C_DOUBLE), pointer :: arr_real(:)
+  complex(C_DOUBLE_COMPLEX), pointer :: arr_complex(:)
+  
+  p_real = fftw_alloc_real(int(npoints, C_SIZE_T))
 
-  p_in = fftw_alloc_real(int(npoints, C_SIZE_T))
+  p_complex = fftw_alloc_complex(int(ncomplex, C_SIZE_T))
 
-  call c_f_pointer(p_in, arr_in, [npoints])
+  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_in(i) = f1(t)
+     arr_real(i) = f1(t)
      
-     write(*,*) t, " ", arr_in(i)
+!     write(*,*) t, " ", arr_real(i)
+     
+  END DO
+
+  call fftw_execute_dft_r2c(plan, arr_real, arr_complex)
+
+  DO i = 1, ncomplex
+
+     write(*,*) abs(arr_complex(i))
      
   END DO
 
-  call fftw_free(p_in)
+  call fftw_free(p_real)
+  call fftw_free(p_complex)
   
 CONTAINS