module mva_resmig_mod
  use  base_struct_mod
  use  base_cstret_mod
  implicit none


  real   , dimension(:), pointer    , private :: kx , ky
  real                              , private :: okz, dkz
  integer                           , private :: nkx, nky, nkz      
  complex, dimension(:), allocatable, private :: tri, tro
  real   , dimension(:), allocatable, private :: str
contains


  integer function mva_resmig_init(kx_in, ky_in, akz, eps) result(stat)
    real, dimension(:), target :: kx_in, ky_in
    type(axis)                 :: akz
    real                       :: eps

    stat = 1
    kx => kx_in
    ky => ky_in

    nkx = size(kx); nky = size(ky)
    nkz = akz%n; okz = akz%o; dkz = akz%d

    stat = base_cstret_init(nkz, nkz, okz, dkz, eps)

write(0,*)"nkx,nky,nkz", nkx, nky, nkz
write(0,*)"okz, dkz", okz, dkz
    allocate(tri(nkz), tro(nkz), str(nkz))

    stat = 0
  end function

  integer function mva_resmig_clos() result(stat)
    stat = 1
    nullify(kx)
    nullify(ky)
    deallocate(tri, tro, str)
    stat = base_cstret_clos()
    stat = 0
  end function

  !!.. rho = s0/s : background slownss / true slowness
  integer function mva_resmig_oper(imgi, imgo, rho, khx, khy) result(stat)
    complex, dimension(:,:,:) :: imgi
    complex, dimension(:,:,:) :: imgo
    real                      :: rho
    real                      :: khx, khy

    integer                   :: ikx, iky, ikz
    real                      :: mks, mkg, kz
    real                      :: kzh, kzm, zzs, zzg

    imgo = 0.

    do iky=1, nky 
    do ikx=1, nkx
       mks = sqrt((kx(ikx)-khx)**2+(ky(iky)-khy)**2)  !!.. |ks|
       mkg = sqrt((kx(ikx)+khx)**2+(ky(iky)+khy)**2)  !!.. |kg|
       do ikz=1, nkz
          kz = okz + (ikz-1)*dkz
          kzh = (4.*kz**2+(mkg-mks)**2)/4.
          kzm = (4.*kz**2+(mkg+mks)**2)/4.
          zzs = 1./rho**2 * kzh*kzm - kz**2 * mks**2
          zzg = 1./rho**2 * kzh*kzm - kz**2 * mkg**2
          if (zzs>=0. .and. zzg>=0) then
             str(ikz) = 0.5/kz*(sqrt(zzs)+sqrt(zzg))
          else
             str(ikz) = -2*dkz
          end if
       end do
       tri = imgi(:,ikx,iky)
       stat = base_cstret_oper(str, tri, tro)
       imgo(:,ikx,iky) = tro
    end do
    end do

    stat = 0
  end function 

end module
