program OW_fftst

use sep
use fftst
use OW_fftu_areal

implicit none

logical                                                     :: fftw
complex,allocatable					:: data(:,:),datam(:,:)
integer						:: iwa,m
integer                                                     :: n1,n2,n3,n4
real                                                        :: d1,d2,d3,d4,o1,o2,o3,o4
complex,allocatable					:: tflds(:)

call sep_init()
call sep_begin_prog()

call from_history("n1",n1);call from_history("o1",o1);call from_history("d1",d1)
call from_history("n2",n2);call from_history("o2",o2);call from_history("d2",d2)
call from_history("n3",n3);call from_history("o3",o3);call from_history("d3",d3)
call from_history("n4",n4);call from_history("o4",o4);call from_history("d4",d4)
write(0,*) n1,n2,n3,n4
call to_history("n1",n1,"out");call to_history("o1",o1,"out");call to_history("d1",d1,"out")
call to_history("n2",n2,"out");call to_history("o2",o2,"out");call to_history("d2",d2,"out")
call to_history("n3",n3,"out");call to_history("o3",o3,"out");call to_history("d3",d3,"out")
call to_history("n4",n4,"out");call to_history("o4",o4,"out");call to_history("d4",d4,"out")

call from_param("fftw",fftw,.true.)
!--------------------------------------------------
! Setting FFTW plans
allocate (tflds(n1));tflds=0
allocate(data(n1,n3));data=0
call owfftws(tflds)
call sreed("in",data,8*n1*n3)
if (fftw) then
   write(0,*) "FFTW"
   do iwa=1,n3
      tflds=data(:,iwa)
      call owfftexes(tflds)
      data(:,iwa)=tflds/sqrt(1.*n1)
   
      tflds=data(:,iwa)
      call owifftexes(tflds)
      data(:,iwa)=tflds/sqrt(1.*n1)
   end do
else
   write(0,*) "FTU"
   m = 1
   do while( m < n1 )
      m = m * 2 
   end do
   allocate(datam(m,n3));datam=0;
   datam(1:n1,:)=data
   call  ft1axis( .true., .true., datam)
   call  ft1axis( .false., .true., datam)
   data=datam(1:n1,:)
end if

call srite("out",data,8*n1*n3)
deallocate(data,tflds)
call sep_end_prog()
end program
