module OW_residmig_mod

use sep
use OW_type
use OW_parms_areal
use OW_cstret

implicit none

contains
!======================================================================
!subroutine ow_resmig(imgi, rho1) 
!
! complex, dimension(:,:,:) :: imgi
! complex :: tri(kz%n), tro(kz%n)
! real    :: str(kz%n)
! real                      :: rho1
! integer                   :: ikx, ish, iz, stat
! real                      :: ks, kg,eps
! real                      :: kz1, kz2, kzs, kzg
!
!   eps=1.
!   stat = base_cstret_init(kz%n, kz%n, kz%o, kz%d, eps)
!
!   do ikx=1,kxm%n
!      do ish=1,ksh%n
!         ks = sqrt((kxm%all(ikx)-ksh%all(ish))**2)
!         kg = sqrt((kxm%all(ikx)+ksh%all(ish))**2)
!         do iz=1,kz%n
!            kz1 = (4.*kz%all(iz)**2+(kg-ks)**2)/4.
!            kz2 = (4.*kz%all(iz)**2+(kg+ks)**2)/4.
!            kzs = 1./rho1**2 * kz1*kz2 - kz%all(iz)**2*ks**2
!            kzg = 1./rho1**2 * kz1*kz2 - kz%all(iz)**2*kg**2
!            if (kzs>=0. .and. kzg>=0) then
!               str(iz) = 0.5/kz%all(iz)*(sqrt(kzs)+sqrt(kzg))
!            else
!               str(iz) = -2*kz%d
!            end if
!         end do
!         tri = imgi(ikx,ish,:)
!         stat = base_cstret_oper(str, tri, tro)
!         imgi(ikx,ish,:) = tro
!      end do
!   end do
!!   stat = base_cstret_clos() 
!end subroutine
!======================================================================
subroutine ow_CAMresmig(imgi, rho1) 

 complex, dimension(:,:,:,:) :: imgi
 complex :: tri(kz%n), tro(kz%n)
 real    :: str(kz%n)
 real                      :: rho1
 integer                   :: ikx, ishx, ishy, iz, stat,iky,jz
 real                      :: ks, kg,eps
 real                      :: kz1, kz2, kzs, kzg, kzx

   eps=1.
   stat = base_cstret_init(kz%n, kz%n, kz%o, kz%d, eps)
   jz=kz%n/2
   do iky=1,kym%n
      do ikx=1,kxm%n
         do ishx=1,kshx%n
            do ishy=1,kshy%n
               ks = sqrt((kxm%all(ikx)-kshx%all(ishx))**2+(kym%all(iky)-kshy%all(ishy))**2)
               kg = sqrt((kxm%all(ikx)+kshx%all(ishx))**2+(kym%all(iky)+kshy%all(ishy))**2)
               do iz=1,kz%n
                  kz1 = (4*kz%all(iz)**2 + (kg-ks)**2)/4
                  kz2 = (4*kz%all(iz)**2 + (kg+ks)**2)/4  

                  kzs = 1./rho1**2 * kz1*kz2 - kz%all(iz)**2 * ks**2
                  kzg = 1./rho1**2 * kz1*kz2 - kz%all(iz)**2 * kg**2

                  if (kzs>=0. .and. kzg>=0) then
                     str(iz) = 0.5/kz%all(iz)*(sqrt(kzs)+sqrt(kzg))
                  else
                     str(iz) = -2*kz%d
                  end if
               end do
               tri = imgi(ikx,ishx,:,iky)
               stat = base_cstret_oper(str, tri, tro)
               imgi(ikx,ishx,:,iky) = tro
            end do

!            do iz=1,kz%n
!               kz1=(kz%all(iz)**2+kym%all(iky)**2+kshx%all(ishx)**2)*&
!                   (kz%all(iz)**2+kym%all(iky)**2+kxm%all(ikx)**2)/&
!                   (kz%all(iz)**2+kym%all(iky)**2)
!               kzs = 1./rho1**2 * kz1 - ks**2
!               kzg = 1./rho1**2 * kz1 - kg**2
!               if (kzs>=0. .and. kzg>=0) then
!                  kzx = 0.5*(sqrt(kzs)+sqrt(kzg))
!                  str(iz) = sqrt(kzx**2 - kym%all(iky)**2)
!               else
!                  str(iz) = -2*kz%d
!               end if
!            end do



         end do
      end do
   end do

end subroutine
end module
