  function z = LinInterp2D(x,y,a,b,n,c,d,m,fA)
%
% Pre:
%    n,m:      integers >=2
%    x,a,b:    scalars that satisfy a<=x<=b 
%    y,c,d:    scalars that satisfy c<=y<=d
%    fA:       an n-by-m matrix with the property that
%              fA(i,j) = f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1))
%              for some function f(.,.). 
% Post:
%    z:        linearly interpolated value for f(x,y).
%   
  hx = (b-a)/(n-1);
  i  = max([1 ceil((x-a)/hx)]);
  dx = (x - (a+(i-1)*hx))/hx;
% x  = a+(i-1+dx)*hx    0<=dx<=1

  hy = (d-c)/(m-1);
  j  = max([1 ceil((y-c)/hy)]); 
  dy = (y - (c+(j-1)*hy))/hy;
% y  = c+(j-1+dy)*hy    0<=dy<=1
   
  z = (1-dy)*((1-dx)*fA(i,j)+dx*fA(i+1,j)) + dy*((1-dx)*fA(i,j+1)+dx*fA(i+1,j+1));
   
   
