program ForFT  
  use sep
  use flags
  use cos_taper_mod 
  implicit none
  integer :: j1,i1,j2,i2,i3,j3,a,b,c,d,i,j,k1,k2,k3,ompi,ompn
  integer :: n1,n2,n3,n4,n5,esize,n1out,blck,nb,er,nf
  real    :: o1,o2,o3,d1,d2,d3,maxf,minf,max1,min1,df,mem,scale
  real   , pointer,dimension(:,:)  :: data
  complex, pointer,dimension(:,:)  :: cdata  
  integer*8,pointer:: planF(:),planI(:)
  logical :: freq
  call sep_init()
  call from_history('n1',n1);call from_history ('n2',n2);call from_history('d1',d1)
  df=1./(d1*n1)
  call from_param('freq',freq,.false.)
  call from_param('maxf',maxf,.5/d1);  call from_param('minf',minf,0.)
  scale=0.075*(maxf-minf)
  call from_param('max1',max1,maxf-scale); call from_param('min1',min1,minf+scale)
  if(freq) call to_history('esize',8)
  call to_history('n1',n1)
  call sep_close()
  nf=n1/2
  allocate(cdata(n1,n2),data(n1,n2))

  call sep_read(data)

write(0,*) 'cdata before plan',sum(cabs(cdata))
write(0,*) ''

  allocate (planF(n2),planI(n2))
write(0,*) 'size(planF),size(planI)',size(planF),size(planI)
  do i=1,n2
     call sfftw_plan_dft_1d(planF(i),n1,&
     cdata(1,i),cdata(1,i),FFTW_FORWARD,FFTW_PATIENT)        
  end do
write(0,*) 'cdata after plan',sum(cabs(cdata))
  if(.not.freq) then
     do i=1,n2
        call sfftw_plan_dft_1d(planI(i),n1,&
        cdata(1,i),cdata(1,i),FFTW_BACKWARD,FFTW_PATIENT)        
     end do
  end if
write(0,*) 'cdata after inverse plan',sum(cabs(cdata))
 write(0,*) 'data',sum(data)  
  do k2=1,n2
     do k1=1,n1
        cdata(k1,k2)=cmplx(data(k1,k2),0.)
     end do
  end do
write(0,*) 'cdata after read',sum(cabs(cdata))
 do k2=1,n2
write(0,*) 'forward'
     call sfftw_execute(planF(k2))
  end do
write(0,*) 'df,dt',df,d1
write(0,*) 'nf',nf
  call init_cos_taper(nf,0.,df,minf,min1,maxf,max1)
  write(0,*) sum(cos_taper_filter)
  do i2=1,n2
     do i1=1,n1
        cdata(i1,i2)=cos_taper_filter(i1)*cdata(i1,i2)
     end do
  end do

  if (.not.freq) then
     do i=1,n2
write(0,*) 'backward'
        call sfftw_execute(planI(i))
     end do
     do k2=1,n2
        do k1=1,n1
           data(k1,k2)=real(cdata(k1,k2))
        end do
     end do
     call sep_write(data)
  else
     call sep_write(cdata)
  end if

 
  write(0,*) 'cdata',sum(cabs(cdata))

end program ForFT
