Program SingleShot
  use sep
  use smooth_mod
!  use cos_taper_mod 
  use flags
  use timers

  implicit none
  integer :: j1,i1,j2,i2,i3,j3,a,b,c,d,i,j,k,k1,k2,k3,ompi,ompn,nof,go,sn
  integer :: n1,n2,n3,n4,n5,esize,n1out,blck,nb,er,nf,nt,trace,rect2,ip
  real    :: o1,o2,o3,d1,d2,d3,df,maxf,minf,mem,scale,min1,max1,dt,small
  real   , pointer,dimension(:,:)  :: data,denom,shotcab
  complex, allocatable,dimension(:,:)  :: cdata,shottrace
  logical:: add,signum,scaled,tapered,decon,one
  integer*8,pointer:: planF(:),planI(:)
  integer, external :: omp_get_thread_num,omp_get_num_threads,sreed_window,srite_window
  real:: startR,finishR,startW,finishW,startFFT,finishFFT, startMult, finishMult,start,finish  
  startR=0.;finishR=0.;startW=0.;finishW=0.;start=0.;finish=0.
  startFFT=0.;finishFFT=0.; startMult=0.; finishMult=0.

  call start_timer(start)
  !----------------------------------------------------------------
  !!     header/parameter section
  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call auxin("traces")
  call from_param('add',add,.false.); call from_param("signum",signum,.false.)

  call from_aux("traces",'n1',n1);call from_aux("traces",'n2',n2);call from_aux("traces",'n3',n3)
  call from_aux("traces",'o1',o1);call from_aux("traces",'o2',o2);call from_aux("traces",'o3',o3)
  call from_aux("traces",'d1',d1);call from_aux("traces",'d2',d2);call from_aux("traces",'d3',d3)
  call from_aux("traces",'esize',esize); if(esize==8) call seperr("need time input")
  call from_param("decon",decon,.false.);call from_param("rect",rect2,1)
  df=1./(d1*n1)
  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)
  call from_param('onlyone',one,.true.)
  call from_param('small',small,0.0001);call from_param("shot",sn,119)
  go=n3
  if(one) go=1
  ompn=1
#ifdef SEP_OMP  
!$OMP PARALLEL 
  ompn = omp_get_num_threads()  
!$OMP END PARALLEL 
#endif
  write(0,*) 'OMP threads',ompn

!  call init_cos_taper(n1/2+1,0.,df,minf,min1,maxf,max1,ompn)
!  bf=i0_n
!  ef=i0_x
  if(add) then
     call auxinout("gathers")
     call from_aux("gathers","esize",er)
     if(er.ne.4) call seperr("can not add frequency and time")
  else 
     call to_history('decon',decon,"gathers");call to_history('rect',rect2,"gathers")
     call to_history('esize',esize,"gathers");call to_history('d1',d1,"gathers")
     call to_history('n1',n1,"gathers");call to_history('o1',o1,"gathers")
     call to_history('n2',n2,"gathers");call to_history('o2',0.,"gathers");call to_history('d2',1.,"gathers")
     call to_history('n3',go,"gathers");call to_history('o3',0.,"gathers");call to_history('d3',1.,"gathers")
     call to_history('n4',1,"gathers");call to_history('o4',0.,"gathers");call to_history('d4',1.,"gathers")
     call to_history('n5',1,"gathers");call to_history('o5',0.,"gathers");call to_history('d5',1.,"gathers")
     call to_history('maxf',minf,"gathers");  call to_history('maxf',maxf,"gathers")
     call to_history('maxf',min1,"gathers");  call to_history('maxf',max1,"gathers")
     call to_history('nt',n1,"gathers")
     call sep_close()
  end if

  write(0,*) 'Sign bit processing:',signum,'   Add :',add
  write(0,*)' n1,n2,n3',n1,n2,n3

  !----------------------------------------------------------------
  !!     allocations
  allocate(data(n1,n2))
  call sep_begin_prog()
  !  mem=mem+4*size(data)
  !  write(0,*) 'allocate data',mem/1024/1024
  allocate(cdata(n1,n2))
  !  mem=mem+8*size(cdata)
  !  write(0,*) 'allocate cdata',mem/1024/1024
  allocate (planF(n2))
  !  mem=mem+8*size(planF)
  !  write(0,*) 'allocate planF',mem/1024/1024
  allocate(shottrace(n1,ompn),shotcab(n1,ompn))
  !  mem=mem+8*size(shottrace)
  !  write(0,*) 'allocate shottrace',mem/1024/1024
  shottrace=(0.,0.)
  allocate (planI(n2))
  !     mem=mem+8*size(planI)
  !     write(0,*) 'allocate planI',mem/1024/1024
  !  write(0,*) 'size(res,3),size(time,1),size(time,2),size(time,3)',size(res,3),size(time,1),size(time,2),size(time,3)
  !----------------------------------------------------------------
  !!    FFT Section
  call start_timer(startFFT)
  do i=1,n2
     call sfftw_plan_dft_1d(planF(i),n1,&
     cdata(1,i),cdata(1,i),FFTW_FORWARD,FFTW_PATIENT)        
  end do
  do i=1,n2
     call sfftw_plan_dft_1d(planI(i),n1,&
     cdata(1,i),cdata(1,i),FFTW_BACKWARD,FFTW_PATIENT)
  end do
  call end_timer(finishFFT,startFFT)
  call sfftw_cleanup()
  allocate(denom(n1,n2))
  !----------------------------------------------------------------  
  events:do ip=1,n3
     write(0,*)' 3 axis: ',ip

     call sep_read(data,"traces")

     if (signum) data=data/abs(data)
     !  call auxclose("traces")
     scale=1./real(n1)
     cdata=cmplx(0.,0.)
     call start_timer(startFFT)

!$OMP PARALLEL DO PRIVATE(i1,i2)
     do i2=1,n2
        do i1=1,n1
           cdata(i1,i2)=cmplx(data(i1,i2),0.)
        end do
        call sfftw_execute(planF(i2))
     end do
!$OMP END PARALLEL DO 
     call end_timer(finishFFT,startFFT)

     do ompi=1,ompn
        do i1=1,n1
           shottrace(i1,ompi)=conjg(cdata(i1,sn))
           shotcab(i1,ompi)=cabs(cdata(i1,sn))
        end do
     end do

     if(decon) then
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
     ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE( i2,i1)
        do i2=1,n2
           do i1=1,n1
              denom(i1,i2)=cabs(cdata(i1,i2))*shotcab(i1,ompi)
              if(denom(i1,i2).lt.small) denom(i1,i2)=1.
!              denom(i1,i2)=cos_taper_filter(i1,ompi)/denom(i1,i2)
              denom(i1,i2)=scale/denom(i1,i2)
           end do
        end do
!$OMP END DO 
!$OMP END PARALLEL
        if(rect2.gt.1)call smooth(denom,1,rect2,0,1,0,0,0)
     else

!write(0,*)' costap',cos_taper_filter
!!$OMP PARALLEL DO PRIVATE(i2,i1)
        do i2=1,n2
           do i1=1,n1
              denom(i1,i2)=scale!*cos_taper_filter(i1,ompi)
           end do
        end do
!!$OMP END PARALLEL DO 
     end if

     ompi=1
!$OMP PARALLEL PRIVATE(ompi)
#ifdef SEP_OMP
     ompi = omp_get_thread_num()+1    
#endif
!$OMP DO PRIVATE( i2,i1)
     do i2=1,n2
        do i1=1,n1
           cdata(i1,i2)=cdata(i1,i2)*denom(i1,i2)*shottrace(i1,ompi)
        end do
        call  sfftw_execute(planI(i2))
     end do
!$OMP END DO 
!$OMP END PARALLEL  

     if(.not.one ) then
        if(add) then
           call sep_read(data,"gathers")
           call sseek("gathers",4*n1*n2*(ip-1),0)
        else
           data=0.
        end if
        do i2=1,n2
           do i1=1,n1
              data(i1,i2)=data(i1,i2)+real(cdata(i1,i2))
           end do
        end do
        
        call sep_write(data,"gathers")
     end if
  end do events
  if(one) then
        if(add) then
           call sep_read(data,"gathers")
           call sseek("gathers",4*n1*n2*(ip-1),0)
        else
           data=0.
        end if
        do i2=1,n2
           do i1=1,n1
              data(i1,i2)=data(i1,i2)+real(cdata(i1,i2))
           end do
        end do
        
        call sep_write(data,"gathers")
     end if

end Program SingleShot
