module OW_scttr_areal

use OW_type
use OW_parms_areal
use OW_fftw_areal

implicit none

contains
!
! From Yaxun library
!
!subroutine owscctr(iwa,sign,ddws,wv,vel)
!
!complex		:: ddws(:), wv(:)
!real		:: vel(:)
!real                :: dx2, w2, s2
!integer             :: sign,ix,iwa
!complex             :: iwdz, lapxo
!
!iwdz = cmplx(0., 1.)*rec%w%all(iwa)*img%z%d*sign 
!dx2  = img%xm%d**2
!w2   = rec%w%all(iwa)**2
!
!do ix=1, kxm%n
!   s2 = vel(ix)
!   s2=s2*s2
!   if (ix == 1) then
!      lapxo = (wv(ix+1)-wv(ix))/img%xm%d
!   else if (ix == kxm%n) then
!      lapxo = (wv(ix)-wv(ix-1))/img%xm%d
!   else
!      lapxo = (wv(ix-1)-2.*wv(ix)+wv(ix+1))/dx2
!   end if
!   ddws(ix) = iwdz*(wv(ix)+0.5/(w2*s2)*lapxo) 
!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(5)
complex			:: tfld(:)
 c=(/1.,0.5,0.375,0.3125,0.2734375/)
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,4
   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*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))
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
