   function u = sht_down_mod(x,z,t,v,wav);
%  Extrapolacao de onda acustica em meio 2D
%  por diferencas finitas 4a ordem em z,x
%
%   function u = we_df_O4(x,z,t,v,wav);
%
% 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);
nwav=length(wav);
dx=abs(x(2)-x(1));
dz=abs(z(2)-z(1));
dt=abs(t(2)-t(1));
kx2=(dt/dx)^2/12;
kz2=(dt/dz)^2/12;
v2=v.^2;

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);

%
% condicoes de contorno e iniciais
%
nxv=round((length(v(1,:))+1)/2);
nxm=round((nx+1)/2);
nxv-nxm+3
nxv-nxm+nx-2
u=zeros(nz,nx,nt-1);
u(1,nxm,2:nwav)=diff(wav);
%
% Inicio dos loops
%
for j=1:nt-1;
    for k=3:nx-2;
        for i=3:nz-2;
%            nxv-nxm+k
            alfax=v(i,nxv-nxm+k)*dt/dx;
            alfaz=v(i,nxv-nxm+k)*dt/dz;
            if k<=2;
                u(i,k,j+1)=alfaz*u(i+1,k,j)+alfax*u(i,k+1,j)+(1-alfax-alfaz)*u(i,k,j);
            elseif k>=(nx-1);
                u(i,k,j+1)=alfaz*u(i+1,k,j)+alfax*u(i,k-1,j)+(1-alfax-alfaz)*u(i,k,j);
            elseif k>2 && k<(nx-1);
                if j==1;
                    X=kz2*v2(i,k)*(-u(i,k+2,j)-u(i,k-2,j)+16*u(i,k+1,j)+16*u(i,k-1,j)-30*u(i,k,j));
                    Z=kx2*v2(i,k)*(-u(i-2,k,j)+16*u(i-1,k,j)+16*u(i+1,k,j)-u(i+2,k,j)-30*u(i,k,j));
                    T=2*u(i,k,j);
                    u(i,k,j+1)=Z+X+T;
                else;
                    X=kz2*v2(i,k)*(-u(i,k+2,j)-u(i,k-2,j)+16*u(i,k+1,j)+16*u(i,k-1,j)-30*u(i,k,j));
                    Z=kx2*v2(i,k)*(-u(i-2,k,j)+16*u(i-1,k,j)+16*u(i+1,k,j)-u(i+2,k,j)-30*u(i,k,j));
                    T=2*u(i,k,j)-u(i,k,j-1);
                    u(i,k,j+1)=Z+X+T;
                end;
            end;
        end;
    end;
%    disp(['tempo = ',num2str(t(j)),' s']);
end;

%for j=1:nt;
%    u(:,:,j)=real(fftshift(ifft2(filt.*fft2(u(:,:,j)))));
%end;
    
end;

