function [ux,uy] = rec_2d_int ...
 ...
  (ax,bx ... % x end points
  ,ay,by ... % y end points
  ,Nx,Ny ... % grid size
  ,glx,gly ... % grid lines
  ,gux,guy ... % grid velocity
  ,x,y ...     % interpolation point
  )

%------------------------------
% Bilinear interpolation of the velocity
% at the point (x, y)
%-------------------------------

%-----------------------------
% locate the x and y intervals
%-----------------------------

   for kx=1:Nx
     prod = (x-glx(kx))*(x-glx(kx+1));
     if(prod<0) break; end
   end

   for ky=1:Ny
     prod = (y-gly(ky))*(y-gly(ky+1));
     if(prod<0) break; end
   end

%----------------------
% interpolation weights
%----------------------

 A00 = (glx(kx+1)-x    )*(gly(ky+1)  -y);
 A01 = (glx(kx+1)-x    )*(y        -gly(ky));
 A10 = (x      -glx(kx))*(gly(ky+1)  -y);
 A11 = (x      -glx(kx))*(y        -gly(ky));

 A = (glx(kx+1)-glx(kx))*(gly(ky+1)-gly(ky));

 w00 = A00/A;
 w01 = A01/A;
 w10 = A10/A;
 w11 = A11/A;

%--------------
% interpolation
%--------------

   ux = w00 * gux(kx  ,ky) ...
      + w01 * gux(kx  ,ky+1) ...
      + w10 * gux(kx+1,ky) ...
      + w11 * gux(kx+1,ky+1);

   uy = w00 * guy(kx  ,ky) ...
      + w01 * guy(kx  ,ky+1) ...
      + w10 * guy(kx+1,ky) ...
      + w11 * guy(kx+1,ky+1);

%---
% done
%----

return