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_resmig_tau(imgi, rho1) 

 complex, dimension(:,:,:) :: imgi
 complex, dimension(:), allocatable :: tri, tro
 real   , dimension(:), allocatable :: str
 real                      :: rho1
 integer                   :: ikx, ish, iz, stat
 real                      :: ks, kg, kz0, kth0,eps
 real                      :: kz1, kz2, kzs, kzg,skt

   allocate(tri(kz%n), tro(kz%n), str(kz%n))
   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
         kth0 = (rho1*s0*ksh%all(ish))**2
         do iz=1,kz%n
            skt=(s0*ksh%all(ish))**2-kz%all(iz)**2
            if (skt >= 0) then
               ks = (kxm%all(ikx)+sqrt(skt))**2
               kg = (kxm%all(ikx)-sqrt(skt))**2
               kz0 = (kz%all(iz)**2+kxm%all(ikx)**2)/kz%all(iz)**2
               kzs = kth0*kz0 - ks
               kzg = kth0*kz0 - kg
            end if
            if (kzs>=0. .and. kzg>=0) then
               str(iz) = 0.5/(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
         if (mod(ikx,100).eq.1.and.mod(ish,10).eq.1) write(0,*) ikx," from ",kxm%n
      end do
   end do
deallocate(tri,tro,str)
end subroutine


!         ks = kxm%all(ikx)**2+(rho1*s0*ksh%all(ish))**2
!         kg = kxm%all(ikx)**2-(rho1*s0*ksh%all(ish))**2
!         kth0 = (rho1*s0*ksh%all(ish))**2
!!         ks = kxm%all(ikx)**2+(2*rho1*s0*ksh%all(ish))**2
!!         kg = kxm%all(ikx)**2-(2*rho1*s0*ksh%all(ish))**2
!!         kth0 = (2*rho1**2*s0*ksh%all(ish))**2
!         do iz=1,kz%n
!            kz0 = (kz%all(iz)**2+kxm%all(ikx)**2)/kz%all(iz)
!            kz1 = ks - kz%all(iz)**2
!            kz2 = kg + kz%all(iz)**2
!            kzs = kth0*kz0 - kz1
!            kzg = kth0*kz0 - kz2
!            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



!======================================================================
subroutine ow_resmig_lin(imgi, rho1) 

 complex, dimension(:,:,:) :: imgi
 real                      :: rho1,rho2
 integer                   :: ikx, ish, iz, stat
 real                      :: ks, kg,mu,ks2,kg2,mu2,eps
 real                      :: kz1, kz2, kzs, kzg,dkzdrho

   eps=1.

   rho2=rho1*rho1

   do ikx=1,kxm%n
      do ish=1,ksh%n
         ks = sqrt((kxm%all(ikx)-ksh%all(ish))**2);ks2=ks*ks
         kg = sqrt((kxm%all(ikx)+ksh%all(ish))**2);kg2=kg*kg
         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.
            if (iz.eq.1) then
               mu=huge(1.);mu2=huge(1.)
               kzs = huge(1.)
               kzg = huge(1.)
               dkzdrho=0.
            else
               mu=kz1*kz2/(16*kz%all(iz));mu2=mu*mu
               kzs = 0.5*sqrt(rho2*mu2-ks2)
               kzg = 0.5*sqrt(rho2*mu2-kg2)
               if ((mu2-ks2.ne.0.).or.(mu2-kg2.ne.0.)) then
                  dkzdrho=0.5*mu2*(1/sqrt(mu2-ks2)+1/sqrt(mu2-kg2))
               end if
            end if
!            write(0,*) ikx,ish,iz,dkzdrho
            imgi(iz,ish,ikx)=imgi(iz,ish,ikx)*dkzdrho
         end do
      end do
   end do

end subroutine
end module
