program FTonesidet2f
  !! 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,nmf,er,ompn,ompi,i
  real::df,maxf,minf,min1,max1,ot,dt,d1,small
  integer*8,pointer:: planF(:)
  complex,  pointer:: cdata(:,:),res(:,:)
  real,     pointer:: data(:,:)
  real::startFFT,finishFFT,startW,finishW,scale,startDIV,finishDIV
  logical::scaleme
  integer, external :: omp_get_thread_num,omp_get_num_threads,sreed_window,srite_window
  call sep_init(SOURCE)
  call from_param('scale',scaleme,.true.)
  call from_param('eps',small,0.001)
  call from_history("n1",nt)
  call from_history("d1",dt)
  call from_history("n2",n2)
  call from_history("n3",n3)
  call from_history("n4",n4)
  call from_history("n5",n5)

!  nt=(nf-1)*2
  nf=nt/2-1
  ot=0.
!  dt=1./(nt*df)
  df=1./(nt*dt)

  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)


  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)
  nmf=i0_x-i0_n+1  
  minf=(i0_n-1)*df
  maxf=(i0_x-1)*df
  if(mod(nmf,2)/=0) nmf=nmf-1

  write(0,'(4x,a,f5.3,a,f6.3)') 'Minimum freq: ',minf,'.  Maximum freq ',maxf
  call to_history("n1",nmf)
  call to_history("d1",df)
  call to_history("o1",minf)
  call to_history("esize",8)
  call to_history("scale",scaleme)
  call to_history("eps",small)
  call sep_close()

  allocate(cdata(nf,n2),data(nt,n2),res(nmf,n2),planF(n2))

  do i2=1,n2
     call sfftw_plan_dft_r2c_1d(planF(i2),nt,&
     data(1,i2),cdata(1,i2),FFTW_PATIENT)        
  end do

  if (.not.scaleme) then
     scale=1./sqrt(real(nt))
  else
     scale=1./real(nt)
  end if
  ompi=1

  sy:do i5=1,n5
     sx:do i4=1,n4
        gy:do i3=1,n3
           call sep_read(data)

           call start_timer(startFFT)
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
              ompi = omp_get_thread_num()+1    
#endif

!$OMP DO PRIVATE(i1,i2,i)
           gx:do i2=1,n2
              call sfftw_execute(planF(i2))
              i1=i0_n
              do i=1,nmf
                 res(i,i2)=cdata(i1,i2)*cos_taper_filter(i1,ompi)*scale
                 i1=i1+1
              end do
           end do gx
!$OMP END DO 
!$OMP END PARALLEL 
           call end_timer(finishFFT,startFFT)           
           call start_timer(startDIV)
!$OMP PARALLEL DO PRIVATE( i1,i2)
           do i2=1,n2
              do i1=1,nmf
                 data(i1,i2)=cabs(res(i1,i2))
                 if(data(i1,i2).gt.small) res(i1,i2)=res(i1,i2)/data(i1,i2)
              end do
           end do
!$OMP END PARALLEL DO 
           call end_timer(finishDIV,startDIV)
           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 Dividing      :', &
  finishDIV,' (s)'
write(0,'(3x,a,f10.3,a)')' Time Writing       :', &
  finishW,' (s)'
deallocate(planF,cdata,data,res)
write(0,*) ' Done'
end program FTonesidet2f
