program FToneside
  !! half the complex axis-- all the time you got comin!
  use sep
  use flags
  use timers  
  use cos_taper_mod
  implicit none

  integer :: i1,i2,i3,i4,i5,n2,n3,n4,n5,nt,nf,nmt,er,ompn,ompi
  real::df,maxf,minf,min1,max1,ot,dt,d1,maxt
  integer*8,pointer:: planI(:)
  complex,  pointer:: cdata(:,:)
  real,     pointer:: data(:,:),res(:,:)
  real::startFFT,finishFFT,startW,finishW,scale
  logical::scaled,tapered
  integer, external :: omp_get_thread_num,omp_get_num_threads,sreed_window,srite_window
  call sep_init(SOURCE)
  call from_param('scaled',scaled,.true.)
  call from_param('tapered',tapered,.true.)

  call from_history("n1",nf)
  call from_history("d1",df)
  call from_history("n2",n2)
  call from_history("n3",n3)
  call from_history("n4",n4)
  call from_history("n5",n5)

  call from_param('maxf',maxf,(nf-1)*df); 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)
  nt=(nf-1)*2
  ot=0.
  dt=1./(nt*df)
  call from_param('maxt',maxt,(nt-1)*dt/2)
  nmt=1.5+maxt/dt
  write(0,'(4x,a,f7.3)') 'Maximum time',maxt
  allocate(cdata(nf,n2),data(nt,n2),res(nmt,n2),planI(n2))
  if (.not.tapered) then
     ompn=1
#ifdef SEP_OMP  
!$OMP PARALLEL 
     ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
#endif
  write(0,'(4x,a,i2 )') 'OMP threads',ompn
     call init_cos_taper(nf,0.,df,minf,min1,maxf,max1,ompn)
  end if
  call to_history("n1",nmt)
  call to_history("d1",dt)
  call to_history("o1",ot)
  call to_history("esize",4)

  do i2=1,n2
     call sfftw_plan_dft_c2r_1d(planI(i2),nt,&
     cdata(1,i2),data(1,i2),FFTW_PATIENT)        
  end do

  if (scaled) then
     scale=1./sqrt(real(nt))
  else
     scale=1./real(nt)
  end if
  sy:do i5=1,n5
     sx:do i4=1,n4
        gy:do i3=1,n3
           call sep_read(cdata)
           if (.not.tapered) then
              ompi=1
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
              ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE(i1,i2)
              do i2=1,n2
                 do i1=1,nf
                    cdata(i1,i2)=cdata(i1,i2)*cos_taper_filter(i1,ompi)
                 end do
              end do
!$OMP END DO 
!$OMP END PARALLEL
           endif
              call start_timer(startFFT)
!$OMP PARALLEL DO PRIVATE( i2,i1)
           gx:do i2=1,n2
              call sfftw_execute(planI(i2))

              w2:do i1=1,nmt
                 res(i1,i2)=data(i1,i2)*scale
              end do w2
           end do gx
!$OMP END PARALLEL DO 
           call end_timer(finishFFT,startFFT)
           call start_timer(startW)
           call sep_write(res)
           call end_timer(finishW,startW)
           er=(i5-1)*n4+(i4-1)*n3+(i3)*n2           
           if(mod(er,(n5*n4*n3*n2)/5)==0) write(0,'(2x,a,i5,a,i5,a,i3,a)') &
             "   Finished",i3,' of',n3*n4*n5,' panels:',100*er/(n5*n4*n3*n2),'%'
        end do gy

     end do sx
  end do sy

write(0,'(3x,a,f10.3,a)')' Time FFTing        :', &
  finishFFT,' (s)'
write(0,'(3x,a,f10.3,a)')' Time Writing       :', &
  finishW,' (s)'
deallocate(planI,cdata,data)
write(0,*) ' Done'
end program FToneside
