   function u = rec_down(x,z,t,v,sism);
%  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

%
% calculos iniciais
%
nx=length(x);
nz=length(z);
nt=length(t);
dx=abs(x(2)-x(1));
dz=abs(z(2)-z(1));
dt=abs(t(2)-t(1));
%
% condicoes de contorno e iniciais
%
u=zeros(nz,nx,nt);
sismd=zeros(size(sism));
sismd(2:nt,:)=diff(sism(:,:));
%
% Inicio dos loops
%
for j=nt:-1:2;
    for i=1:nx;
        u(1,i,j)=u(1,i,j)+sismd(j,i);
        for k=1:nz;
            if j==nt;
                T=u(k,i,j);
            else;
                T=u(k,i,j+1)-2*u(k,i,j);
            end;
            
            if i==1;
                X=dx^(-2)*(u(k,i+1,j)-2*u(k,i,j));
            elseif i==nx;
                X=dx^(-2)*(-2*u(k,i,j)+u(k,i-1,j));
            else;
                X=dx^(-2)*(u(k,i+1,j)-2*u(k,i,j)+u(k,i-1,j));
            end;
            
            if k==1;
                Z=dz^(-2)*(u(k+1,i,j)-2*u(k,i,j));
            elseif k==nz;
                Z=dz^(-2)*(-2*u(k,i,j)+u(k-1,i,j));
            else;
                Z=dz^(-2)*(u(k+1,i,j)-2*u(k,i,j)+u(k-1,i,j));
            end;
            
            u(k,i,j-1)=-T+(Z+X)*(v(k,i)*dt)^(2);
            
        end;
    end;
%    disp(['tempo = ',num2str(t(j)),' s']);
end;
