   function u = rec_down(xs,is,xh,xm,z,t,v,s1,npd);
%  Extrapolacao de onda acustica em meio 2D
%  por diferencas finitas 4a ordem em z,x
%
%   function u = rec_down(x,z,t,v,sism);
%
% wav   vetor   nwav X 1
%  v    matriz    nx X nz
%  x    vetor     nx X 1
%  z    vetor     nz X 1 
%  t    vetor     nt X 1

%
% Preliminary computations
%
ns=length(xs);
nx=length(xh)+2*npd;
nz=length(z);
nt=length(t);

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.004)^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
%
u=zeros(nz,nx,nt);
sismd=zeros(nt,nx);
% for is=1:ns;
%     sismd(2:nt,npd+1:npd+size(s1,2),is)=diff(s1(:,:,is));
% end;
     sismd(2:nt,npd+1:npd+size(s1,2))=diff(s1(:,:));
%
% Indexing velocity
%
ni=floor((xs(is)-npd*dx-abs(min(xh))-oxm)/dx)+1;
vel=v(:,ni:ni+nx-1);
%
% Loops
%
for j=nt:-1:2;
    for k=1:nx;
        u(1,k,j)=u(1,k,j)+sismd(j,k);
        for i=1:nz;
            Z=0;X=0;T=0;
            alfa=vel(i,k)*dt;
            
            if j==nt;
                T=u(i,k,j);
            else;
                T=u(i,k,j+1)-2*u(i,k,j);
            end;
            
            if k>4 && k<(nx-3);
                X=dx^(-2)*(-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==1;
                Z=dz^(-2)*(u(i+1,k,j)-2*u(i,k,j));
            elseif i==2 ;
                Z=dz^(-2)*(u(i+1,k,j)+u(i-1,k,j)-2*u(i,k,j));
            elseif i==3 || i==4;
                Z=dz^(-2)/12*(-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=dz^(-2)*(-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;
            
            u(i,k,j-1)=-T+(Z+X)*(vel(i,k)*dt)^(2);
            
%            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;
    end;
%    disp(['tempo = ',num2str(t(j)),' s']);
%end;
