function [Axx,Axy,Ayx,Ayy] = rec_2d_vgt ...
...
      (glx,gly,Nx,Ny,gux,guy)

%------------------------------
% compute the velocity gradient
% tensor A_ij at the grid points
% by numerical differentiation
%-------------------------------

%----------------
% interior points
% compute derivatives by central differences
%----------------

for i=2:Nx
 for j=2:Ny
   Axx(i,j) = (gux(i+1,j)-gux(i-1,j))/(glx(i+1)-glx(i-1));
   Axy(i,j) = (guy(i+1,j)-guy(i-1,j))/(glx(i+1)-glx(i-1));
   Ayx(i,j) = (gux(i,j+1)-gux(i,j-1))/(gly(j+1)-gly(j-1));
   Ayy(i,j) = (guy(i,j+1)-guy(i,j-1))/(gly(j+1)-gly(j-1));
 end
end

%----------
% left wall:
% compute derivatives by central or forward differences
%----------

i=1;
for j=2:Ny
  Axx(i,j) = (gux(i+1,j)-gux(i,j))/(glx(i+1)-glx(i));
  Axy(i,j) = (guy(i+1,j)-guy(i,j))/(glx(i+1)-glx(i));
  Ayx(i,j) = (gux(i,j+1)-gux(i,j-1))/(gly(j+1)-gly(j-1));
  Ayy(i,j) = (guy(i,j+1)-guy(i,j-1))/(gly(j+1)-gly(j-1));
end

%----------
% bottom wall:
% compute derivatives by central or forward differences
%----------

j=1;
for i=2:Nx
  Axx(i,j) = (gux(i+1,j)-gux(i-1,j))/(glx(i+1)-glx(i-1));
  Axy(i,j) = (guy(i+1,j)-guy(i-1,j))/(glx(i+1)-glx(i-1));
  Ayx(i,j) = (gux(i,j+1)-gux(i,j))/(gly(j+1)-gly(j));
  Ayy(i,j) = (guy(i,j+1)-guy(i,j))/(gly(j+1)-gly(j));
end

%----------
% right wall:
% compute derivatives by central or backward differences
%----------

i=Nx+1;
for j=2:Ny
 Axx(i,j) = (gux(i,j)-gux(i-1,j))/(glx(i)-glx(i-1));
 Axy(i,j) = (guy(i,j)-guy(i-1,j))/(glx(i)-glx(i-1));
 Ayx(i,j) = (gux(i,j+1)-gux(i,j-1))/(gly(j+1)-gly(j-1));
 Ayy(i,j) = (guy(i,j+1)-guy(i,j-1))/(gly(j+1)-gly(j-1));
end

%----------
% top wall:
% compute derivatives by central or backward differences
%----------

j=Ny+1;
for i=2:Nx
 Axx(i,j) = (gux(i+1,j)-gux(i-1,j))/(glx(i+1)-glx(i-1));
 Axy(i,j) = (guy(i+1,j)-guy(i-1,j))/(glx(i+1)-glx(i-1));
 Ayx(i,j) = (gux(i,j)  -gux(i,j-1))/(gly(j)  -gly(j-1));
 Ayy(i,j) = (guy(i,j)  -guy(i,j-1))/(gly(j)  -gly(j-1));
end

%-------------------
% four corner points
%-------------------

i=1;
j=1;
Axx(i,j) = (gux(i+1,j)-gux(i,j))/(glx(i+1)-glx(i));
Axy(i,j) = (guy(i+1,j)-guy(i,j))/(glx(i+1)-glx(i));
Ayx(i,j) = (gux(i,j+1)-gux(i,j))/(gly(j+1)-gly(j));
Ayy(i,j) = (guy(i,j+1)-guy(i,j))/(gly(j+1)-gly(j));
i=Nx+1;
j=1;
Axx(i,j) = (gux(i,j)-gux(i-1,j))/(glx(i)-glx(i-1));
Axy(i,j) = (guy(i,j)-guy(i-1,j))/(glx(i)-glx(i-1));
Ayx(i,j) = (gux(i,j+1)-gux(i,j))/(gly(j+1)-gly(j));
Ayy(i,j) = (guy(i,j+1)-guy(i,j))/(gly(j+1)-gly(j));
i=Nx+1;
j=Ny+1;
Axx(i,j) = (gux(i,j)-gux(i-1,j))/(glx(i)-glx(i-1));
Axy(i,j) = (guy(i,j)-guy(i-1,j))/(glx(i)-glx(i-1));
Ayx(i,j) = (gux(i,j)-gux(i,j-1))/(gly(j)-gly(j-1));
Ayy(i,j) = (guy(i,j)-guy(i,j-1))/(gly(j)-gly(j-1));
i=1;
j=Ny+1;
Axx(i,j) = (gux(i+1,j)-gux(i,j))/(glx(i+1)-glx(i));
Axy(i,j) = (guy(i+1,j)-guy(i,j))/(glx(i+1)-glx(i));
Ayx(i,j) = (gux(i,j)-gux(i,j-1))/(gly(j)-gly(j-1));
Ayy(i,j) = (guy(i,j)-guy(i,j-1))/(gly(j)-gly(j-1));

%-----
% done
%-----

return