 module sht_down_mod

! (xs,is,xh,xm,z,t,v,wav,it0,npd)
!  Extrapolacao de onda acustica em meio 2D
!  por diferencas finitas 8a ordem em z,x
!
!   function u = sht_down_2t8x(xs,is,xh,z,t,v,wav,it0,npd)
!
! wav   vetor     nwav X 1
!  v    matriz    nx X nz
!  xs   vetor     nx X 1        --> shot positions
!  xh   vetor     nx X 1        --> offset vector
!  z    vetor     nz X 1 
!  t    vetor     nt X 1
!  npd  scalar    must be odd   --> # traces to be added on both sides
!  it0  scalar    must be odd 

 implicit none

 contains

subroutine shot_ext(xs,is,xh,xm,z,t,v,wav,it0,npd)

real, allocatable			:: taper(:,:), u(:,:,:), u1(:,:,:)
integer					:: ns,nx,nz,nt,nwav,nxm
integer					:: i,j,k
real					:: dx,dz,dt,oxm,kx2,kz2,
real					:: X,T,Z,alfa

!
! Preliminary computations
!
ns=length(xs)
nx=length(xh)+2*npd
nz=length(z)
nt=length(t)
nwav=length(wav)

dx=abs(xh(2)-xh(1))
dz=abs(z(2)-z(1))
dt=abs(t(2)-t(1))

oxm=min(xm)
kx2=(dt/dx)**2
kz2=(dt/dz)**2

allocate(taper(nz,nx),u(nz,nx,3), u1(nz,nx,3))

taper=1.
u=0.
(nz,nx,it0+nt-1)


!
! Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
!                   SEG Expanded Abstracts 23, 1921 (2004)
!
do i=nz-9:nz
    taper(i,:)=exp(((-1)*(i-nz+9)*0.008)**2)
end do
do i=1:10
    taper(:,i)=exp(((-1)*(10-i)*0.004)**2)
end do
do i=nx-9:nx
    taper(:,i)=exp(((-1)*(i-nx+9)*0.004)**2)
end do
!
! Initial conditions
!
nxm=npd+abs(min(xh))/dx+1
u(1,nxm,1:nwav-1)=diff(wav)
!
! Indexing velocity
!
ni=int((xs(is)-npd*dx-abs(min(xh))-oxm)/dx)+1
vel=v(:,ni:ni+nx-1)
vel2=vel**2
!
! Loops
!
do j=2:it0+nt-1
    do k=1:nx
        do i=2:nz
            Z=0;X=0;T=0
            alfa=vel(i,k)*dt
!
! 8th order propagator in space: Fornberg, B - Mathematics of Computation, Vol. 51, No. 184. (Oct., 1988)
!
            if (k>4 && k<(nx-3)) then
                X=kx2*vel2(i,k)*(-0.00178571428571*(u(i,k+4,j)+u(i,k-4,j))+ & 
                                0.02539682539683*(u(i,k+3,j)+u(i,k-3,j))- &
                                0.20000000000000*(u(i,k+2,j)+u(i,k-2,j))+ &
                                1.60000000000000*(u(i,k+1,j)+u(i,k-1,j))- &
                                2.84722222222222*(u(i,k,j)))
            end if

            if (i=2) then
                Z=kz2*vel2(i,k)*(u(i+1,k,j)+u(i-1,k,j)-2*u(i,k,j))
            else if (i=3 || i=4) then
                Z=kz2/12*vel2(i,k)*(-u(i+2,k,j)-u(i-2,k,j)+16*u(i+1,k,j)+16*u(i-1,k,j)-30*u(i,k,j))
            else if (i>4 && i<(nz-3)) then
                Z=kz2*vel2(i,k)*(-0.00178571428571*(u(i+4,k,j)+u(i-4,k,j))+ & 
                                0.02539682539683*(u(i+3,k,j)+u(i-3,k,j))- &
                                0.20000000000000*(u(i+2,k,j)+u(i-2,k,j))+ &
                                1.60000000000000*(u(i+1,k,j)+u(i-1,k,j))- &
                                2.84722222222222*(u(i,k,j)))
            end if
            
            T=2*u(i,k,j)-u(i,k,j-1)

            u(i,k,j+1)=(Z+X+T)*taper(i,k)

!
! Boundary conditions: Reynolds, A. - GEOPHYSICS. VOL. 43, NO. 6 (OCTOBER 1978)
!
            if (k<=4) then
               u(i,k,j+1)=u(i,k,j)+u(i,k+1,j)-u(i,k+1,j-1)+ &
                   alfa/dx*(u(i,k+1,j)-u(i,k,j)-(u(i,k+2,j-1)-u(i,k+1,j-1)))
               u(i,k,j+1)=taper(i,k)*u(i,k,j+1)
            else if (k>=nx-3) then
               u(i,k,j+1)=u(i,k,j)+u(i,k-1,j)-u(i,k-1,j-1)- &
                   alfa/dx*(u(i,k,j)-u(i,k-1,j)-(u(i,k-1,j-1)-u(i,k-2,j-1)))
               u(i,k,j+1)=taper(i,k)*u(i,k,j+1)
            end if
            
            if (i>=nz-3) then
               u(i,k,j+1)=u(i,k,j)+u(i-1,k,j)-u(i-1,k,j-1)- &
                   alfa/dz*(u(i,k,j)-u(i-1,k,j)-(u(i-1,k,j-1)-u(i-2,k,j-1)))   
               u(i,k,j+1)=taper(i,k)*u(i,k,j+1)
            end if
        end do
    end do
end do
u1=u(:,:,it0:it0-1+nt)
end subroutine
end module
