function [z]=BICUBICINTERPOL(xdata,ydata,zdata,dxdata,dydata,dxdydata,x,y)
%BICUBICINTERPOL function perform a bicubic interpolation along the plane 
% xy for the underlying function z at data point with coordinates x and y.
% The dxdata represents the partial derivatives of the underlying function
% at points on xdata and ydata, dydata represents the partial derivatives
% of the underlying function at points on xdata and ydata, and dxdydata are
% the cross double derivatives at the same points.
%   Data values on xdata, ydata and zdata must surround xy point.
%   BICUBICINTERPOL perform extrapolations as well.
%
%
%                   (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
% dxdata=[dx1,dx2,dx3,dx4]
% dydata=[dy1,dy2,dy3,dy4]
% dzdata=[dz1,dz2,dz3,dz4]
% dxdydata=[dxdy1,dxdy2,dxdy3,dxdy4]
%
%
%

EX=(x-xdata(1))/(xdata(2)-xdata(1));
EY=(y-ydata(2))/(ydata(3)-ydata(2));


F1=((-1+EX)^2)*(1+2*EX)*((-1+EY)^2)*(1+2*EY);
F2=(-EX^2)*(-3+2*EX)*((-1+EY)^2)*(1+2*EY);
F3=(EX^2)*(-3+2*EX)*(EY^2)*(-3+2*EY);
F4=-((-1+EX)^2)*(1+2*EX)*(EY^2)*(-3+2*EY);
F5=((-1+EX)^2)*(EX)*((-1+EY)^2)*(1+2*EY);
F6=(-1+EX)*(EX^2)*((-1+EY)^2)*(1+2*EY);
F7=-(-1+EX)*(EX^2)*(EY^2)*(-3+2*EY);
F8=-((-1+EX)^2)*(EX)*(EY^2)*(-3+2*EY);
F9=((-1+EX)^2)*(1+2*EX)*((-1+EY)^2)*EY;
F10=-(EX^2)*(-3+2*EX)*((-1+EY)^2)*EY;
F11=-(EX^2)*(-3+2*EX)*(-1+EY)*(EY^2);
F12=((-1+EX)^2)*(1+2*EX)*(-1+EY)*(EY^2);
F13=((-1+EX)^2)*(EX)*((-1+EY)^2)*EY;
F14=(-1+EX)*(EX^2)*((-1+EY)^2)*EY;
F15=(-1+EX)*(EX^2)*(-1+EY)*(EY^2);
F16=((-1+EX)^2)*(EX)*(-1+EY)*(EY^2);

z=F1*zdata(1)+F2*zdata(2)+F3*zdata(3)+F4*zdata(4)+...
    F5*dxdata(1)+F6*dxdata(2)+F7*dxdata(3)+F8*dxdata(4)+...
    F9*dydata(1)+F10*dydata(2)+F11*dydata(3)+F12*dydata(4)-...
    F13*dxdydata(1)-F14*dxdydata(2)-F15*dxdydata(3)-F16*dxdydata(4);

end