add second input function for dft

3 files changed, 123 insertions, 88 deletions

diff --git a/Makefile b/Makefile
index 217f0d0..95c1dad 100644
--- a/Makefile
+++ b/Makefile
@@ -2,21 +2,21 @@ FC = gfortran
 FFLAGS = -std=f2008 -Wall -pedantic -fbounds-check -fimplicit-none \
 	 -I/usr/include -lfftw3 -lm
 
-all : f1
+all : fourier
 
-f1 : src/f1.f90
+fourier : src/fourier.f90
 	$(FC) $(FFLAGS) $^ -o $@
 
-res/f1_original.txt : f1
-	./f1 > res/f1_original.txt
+res/f1_original.txt : fourier
+	./fourier f1 > res/f1_original.txt
 
-res/f1_dft.txt : f1
-	./f1 dft > res/f1_dft.txt
+res/f1_dft.txt : fourier
+	./fourier f1 dft > res/f1_dft.txt
 
 res/%.png : src/%.gnuplot res/%.txt
 	gnuplot -c $<
 
 clean :
-	-rm f1 res/*.txt
+	-rm fourier res/*.txt
 
 .PHONY : all clean 
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