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
