function [z]=BILINEARINTERPOL(xdata,ydata,zdata,x,y)
%BILINTERPOL function perform a bilinear interpolation along the plane xy
%for the underlying function z at data point with coordinates x and y.
%   Data values on xdata, ydata and zdata must surround xy point
%
%                   (x4,y4)          (x3,y3)
%                (0,1) * -------------- * (1,1)
%                      |                |
%                      |                |
%                      |        *       |
%                      |      (x,y)     |
%                      |                |
%                (0,0) * -------------- * (1,0)
%                   (x1,y1)          (x2,y2)
%
%  
% xdata=[x1,x2,x3,x4]
% ydata=[y1,y2,y3,y4]
% zdata=[z1,z2,z3,z4]
% x1=x4
% x2=x3
% y1=y2
% y3=y4

r=(x-xdata(1))/(xdata(2)-xdata(1));
s=(y-ydata(3))/(ydata(2)-ydata(3));
z=(1-r)*(1-s)*zdata(4)+r*(1-s)*zdata(3)+(1-r)*s*zdata(1)+r*s*zdata(2);

end