   function u1 = sht_down_2t8x(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 

%
% 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));
if ns~=1;
    ds=abs(xs(2)-xs(1));
else;
    ds=0;
end;

oxm=min(xm);
kx2=(dt/dx)^2;
kz2=(dt/dz)^2;
%
% Post propagation filter
%
nz2=floor(nz/2)+1;
nx2=floor(nx/2)+1;
filt=zeros(nz,nx);
filt(nz2,nx2)=1;   filt(nz2+1,nx2)=1;   filt(nz2+2,nx2)=1;
filt(nz2,nx2+1)=1; filt(nz2+1,nx2+1)=4; filt(nz2+2,nx2+1)=1;
filt(nz2,nx2+2)=1; filt(nz2+1,nx2+2)=1; filt(nz2+2,nx2+2)=1;
filt=fft2(filt);
%
% Taper Sponge-like: Bording, P - Finite difference modeling - nearly optimal sponge boundary conditions
%                   SEG Expanded Abstracts 23, 1921 (2004)
%
taper=ones(nz,nx);
for iz=nz-9:nz;
    taper(iz,:)=exp(((-1)*(iz-nz+9)*0.008)^2);
end;
for ix=1:10;
    taper(:,ix)=exp(((-1)*(10-ix)*0.004)^2);
end;
for ix=nx-9:nx;
    taper(:,ix)=exp(((-1)*(ix-nx+9)*0.004)^2);
end;
%
% Initial conditions
%
nxm=npd+abs(min(xh))/dx+1;
u=zeros(nz,nx,it0+nt-1);
t=((-it0+1)*dt:dt:(it0+nt-1)*dt);
u(1,nxm,1:nwav-1)=diff(wav);
%
% Indexing velocity
%
ni=floor((xs(is)-npd*dx-abs(min(xh))-oxm)/dx)+1;
vel=v(:,ni:ni+nx-1);
vel2=vel.^2;
%
% Loops
%
for j=2:it0+nt-1;
    for k=1:nx;
        for 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);
                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 i==2 ;
                Z=kz2*vel2(i,k)*(u(i+1,k,j)+u(i-1,k,j)-2*u(i,k,j));
            elseif i==3 || i==4;
                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));
            elseif i>4 && i<(nz-3);
                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;
            
            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;
               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);
            elseif k>=nx-3;
               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 i>=nz-3;
               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;
        end;
    end;
    
%    disp(['tempo = ',num2str(t(j)),' s']);
end;
u1=u(:,:,it0:it0-1+nt);
%for j=1:nt;
%    u1(:,:,j)=real(fftshift(ifft2(filt.*fft2(u(:,:,it0-1+j)))));
%end;