module OW_scttr_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal

implicit none

contains
!
! From Yaxun library
!
subroutine owscctrx(ww,sign,ddws,wv,vel)

complex		:: ddws(:,:), wv(:,:), ws(kxm%n+2,kym%n+2)
real		:: vel(:,:)
real                :: dx2, dy2, w2, s2,ww
integer             :: sign,ix,iy,iwa
complex             :: iwdz, lapxo

iwdz = cmplx(0., 1.)*ww*sign*img%z%d
dx2  = img%xm%d**2
dy2  = img%ym%d**2
w2   = ww**2
ws=0;ws(2:kxm%n+1,2:kym%n+1)=wv

do iy=2, kym%n+1
do ix=2, kxm%n+1
   s2 = vel(ix-1,iy-1)
   s2=s2*s2
   lapxo = (ws(ix-1,iy)-2.*ws(ix,iy)+ws(ix+1,iy))/dx2+(ws(ix,iy-1)-2.*ws(ix,iy)+ws(ix,iy+1))/dy2
   ddws(ix-1,iy-1) = iwdz*(ws(ix,iy)+0.5/(w2*s2)*lapxo) 
end do
end do
end subroutine
!
! Paul Sava and Ioan Vlad paper
!!---------------------------------------------------------
!subroutine owscctr(iwa,inode,sign,ddws,wv,vel,tfld,gat)
!character*3		:: gat
!integer			:: sign, j, iwa, inode
!complex			:: dwss(kxm%n),wv(:),ddws(:),dwkk(kxm%n)
!real			:: vel(:),c(10)
!complex			:: tfld(:)
! c=(/1.,1/2,3/8,5/16,35/128,63/256,231/1024,429/2048,6435/32768,12155/65536/)
!if (gat.eq."sou") then
!   tfld=wv
!   call owfftexes(tfld,inode)
!   dwkk=tfld/sqrt(1.*kxm%n)
!else
!   tfld=wv
!   call owfftexer(tfld,inode)
!   dwkk=tfld/sqrt(1.*kxm%n)
!end if
!do j=0,7
!   dwss=sign*dwkk*abs(kxm%all)**(2*j)
!   if (gat.eq."sou") then
!      tfld=dwss
!      call owifftexes(tfld,inode)
!      dwss=tfld/sqrt(1.*kxm%n)*(1/(vel*rec%w%all(iwa)))**(2*j)*c(j+1)
!   else
!      tfld=dwss
!      call owifftexer(tfld,inode)
!      dwss=tfld/sqrt(1.*kxm%n)*(1/(vel*rec%w%all(iwa)))**(2*j)*c(j+1)
!   end if
!   ddws=ddws+dwss
!end do
!ddws=ddws*rec%w%all(iwa)*exp(cmplx(0,pi/2))*img%z%d
!end subroutine
!!---------------------------------------------------------
!subroutine owscctr1(iwa,inode,sign,ddws,wv,vel,tfld,gat)
!character*3		:: gat
!integer			:: sign, j, iwa, inode
!complex			:: dwss(kxm%n),wv(:),ddws(:),dwkk(kxm%n)
!real			:: vel(:),c(10)
!complex			:: tfld(:)
! c=(/1.,1/2,3/8,5/16,35/128,63/256,231/1024,429/2048,6435/32768,12155/65536/)
!if (gat.eq."sou") then
!   tfld=wv
!   call owfftexes(tfld,inode)
!   dwkk=tfld/sqrt(1.*kxm%n)
!else
!   tfld=wv
!   call owfftexer(tfld,inode)
!   dwkk=tfld/sqrt(1.*kxm%n)
!end if
!do j=0,9
!   dwss=sign*dwkk*abs(kxm%all)**(2*j)
!   if (gat.eq."sou") then
!      tfld=dwss
!      call owifftexes(tfld,inode)
!      dwss=tfld/sqrt(1.*kxm%n)*(vel/rec%w%all(iwa))**(2*j)*c(j+1)
!   else
!      tfld=dwss
!      call owifftexer(tfld,inode)
!      dwss=tfld/sqrt(1.*kxm%n)*(vel/rec%w%all(iwa))**(2*j)*c(j+1)
!   end if
!   ddws=ddws+dwss
!end do
!ddws=ddws*rec%w%all(iwa)!*img%z%d
!end subroutine
!
!!---------------------------------------------------------
!!subroutine owscctr(iwa,inode,sign,ddws,wv,vel,tfld,gat)
!!character*3		:: gat
!!integer			:: sign, j, iwa, inode
!!complex			:: dwss(kxm%n),wv(:),ddws(:)
!!real			:: vel(:),c(5)
!!complex			:: tfld(:)
!! c=(/1.,0.5,0.375,0.3125,0.2734375/)
!!do j=0,1
!!   dwss=sign*abs(kxm%all)**(2*j)*wv
!!   if (gat.eq."sou") then
!!      tfld=dwss
!!      call owifftexes(tfld,inode)
!!      dwss=tfld/sqrt(1.*kxm%n)*(vel/rec%w%all(iwa))**(2*j)*c(j+1)
!!   else
!!      tfld=dwss
!!      call owifftexer(tfld,inode)
!!      dwss=tfld/sqrt(1.*kxm%n)*(vel/rec%w%all(iwa))**(2*j)*c(j+1)
!!   end if
!!   ddws=ddws+dwss
!!end do
!!ddws=ddws*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
!!end subroutine
end module
