module RWE2D_WI_wemig
  use sep
  use util_ 
  use fft 
  use Kernel_TOPO

  implicit none
  integer, private :: nt,ng,nw,nr,nzcut,ntap
  real,    private :: dt,dg,dw,ot,og,ow,okg,dkg,eps,pi
  logical, private :: kin,verbose

  complex, allocatable,private :: rax(:,:)
  real,    allocatable,private :: vmin(:),veltmp(:,:),tap(:)
 
contains

  !----------------------------------------------------------------  
  subroutine wemig_init(ag,at,aw,ar,verbose_in,nzcut_in)
    type(myaxis)     :: ag,at,aw,ar
    integer          ::                        nzcut_in
    logical          ::              verbose_in

    integer :: ig
    pi=acos(-1.)
    nt=at%n;    ng=ag%n;    nw=aw%n
    dt=at%d;    dg=ag%d;    dw=aw%d * 2.*pi
    ot=at%o;    og=ag%o;    ow=aw%o * 2.*pi
    nzcut=nzcut_in

    okg=  -pi/    dg
    dkg=2.*pi/(ng*dg)

    verbose=verbose_in;    
    eps=0.0000001

    allocate (tap(ng) )
    allocate(rax(ng,1),vmin(nt),veltmp(ng,1) )
    tap=1.
    call from_param("ntap",ntap,10)
    call rwetaper_init()

  end subroutine wemig_init

  !----------------------------------------------------------------

  !! WE migration driver
  subroutine wemig(rwf,injector,wfld,adj,iss,fields,gnorm,w,RVel)
    integer          :: iw,it,ig,ir,iss
    real             :: rarg,w
    real             :: fields(:,:,:),gnorm(:,:),RVel(:,:)
    complex          :: rwf(:,:),injector(:,:)
    complex          :: wfld(:,:)
    logical          :: logic,adj

    call shot_ker_init(ng, 1, dg , 1. ,dt , 0., 0., .122996)

    !! Get minvel for phase filter
    do it=1,nt
       vmin(it)=minval(RVel(it,:))
    end do

    do ig=1,ng
       rax(ig,1) = rwf(ig,1) 
    end do



    !! Taper at ends
    call rwetaper(rax(:,1)) 

    if (adj) then

       !! Store entire monochromatic wavefield
       do ig=1,ng
          wfld(ig,1)=rax(ig,1)
       end do

       do it=2,nt

          !! SEt up veltmp
          do ig=1,ng
             veltmp(ig,1)=RVel(ig,it)
          end do

          !! Inject wavefield
          rax(nzcut,1) = rax(nzcut,1)+injector(it-1,1)

          !! Propagate one step
          call shot_ker_onestep(rax,veltmp,fields(:,:,it),w,-1.)

          !! Fourier Transform
          call fth(.false.,.false.,rax(:,1) )

          !! Phase correction
          call wemig_phs_correction(rax,w,vmin(it),-1.)

          !! Fourier Transform back
          call fth( .true.,.false.,rax(:,1) )

          !! Taper
          call rwetaper(rax(:,1)) 

          !! Save into total wavefield
          do ig=1,ng
             wfld(ig,it)=rax(ig,1)
          end do
       end do !! . . z-loop

    else !! Propagate in reverse

       !! Store wavefield
       wfld(nzcut,nt)=injector(nt,1)

       do it=nt,2,-1

          !! SEt up veltmp
          do ig=1,ng
             veltmp(ig,1)=RVel(ig,it)
          end do

          !! Inject wavefield
          rax(nzcut,1)=rax(nzcut,1)+injector(it,1)

          !! Propagate one step
          call shot_ker_onestep(rax,RVel,fields(:,:,it),w,1.)

          !! Fourier Transform
          call fth(.false.,.false.,rax(:,1) )

          !! Phase correction
          call wemig_phs_correction(rax,w,vmin(it),1.)

          !! Fourier Transform back
          call fth( .true.,.false.,rax(:,1) )

          !! Taper
          call rwetaper(rax(:,1)) 

          !! Save into total wavefield
          do ig=1,ng
             wfld(ig,it-1)=rax(ig,1) 
          end do

       end do !! . . z-loop
    end if

    call shot_ker_release()

  end subroutine wemig

!----------------------------------------------------------------


  subroutine rwetaper_init()
    integer :: itap,jtap,ix
    real :: pi
    pi=acos(-1.)

    if (ntap .gt. 1) then

       do itap=1,ntap
          jtap = abs(ntap-itap)
          tap(   itap) = cos(pi/2.* jtap/(ntap-1))
       end do

       do itap=1,ntap
          jtap = abs(ntap-itap)
          tap(ng-itap+1) = cos(pi/2.* jtap/(ntap-1))
       end do

    end if

  end subroutine rwetaper_init

!----------------------------------------------------------------

  subroutine rwetaper(dax)
    complex  :: dax(:)
    integer :: ig
    do ig=1,ng
       dax(ig) = dax (ig)* tap(ig)
    end do
  end subroutine rwetaper

end module  RWE2D_WI_wemig
