%
%This function does a "near-bilinear interpolation".
%
%We have to work with distorted quadrilaterals rather than rectangles, but 
%we use the distance to the quadrilateral sides as in bilinear interpolation:
%-> when the quadrilateral is a rectangle this actually is bilinear interpolation.
%-> for convex quadrilaterals it is close.
%-> for non-convex quadrilaterals !? 
%
%Author: Gael Forget
%Date: June 14th 2007
%
function [coeffs]=netcdf_ecco_GenericgridCoeffs(x0,y0,x_all,y_all);

warning off MATLAB:divideByZero;

x=ones([length(x_all) size(x0)]); for icur=1:length(x_all); x(icur,:,:)=x0; end;
y=ones([length(x_all) size(x0)]); for icur=1:length(x_all); y(icur,:,:)=y0; end;

alpha=zeros(size(x)); beta=zeros(size(x)); gamma=zeros(size(x)); coeffs=zeros(size(x));

x_all2=zeros([length(x_all) size(x0,1) 1]); for icur=1:length(x_all); x_all2(icur,:,:)=x_all(icur); end;
y_all2=zeros([length(y_all) size(x0,1) 1]); for icur=1:length(y_all); y_all2(icur,:,:)=y_all(icur); end;

for icur=1:4;

x=circshift(x,[0 0 1]); y=circshift(y,[0 0 1]);
alpha=circshift(alpha,[0 0 1]); beta=circshift(beta,[0 0 1]); 

v1x=x_all2(:,:)-x(:,:,1); v1y=y_all2(:,:,1)-y(:,:,1); v2x=x(:,:,2)-x(:,:,1); v2y=y(:,:,2)-y(:,:,1);
alpha(:,:,1)=(v1x.*v2x+v1y.*v2y)./(v2x.*v2x+v2y.*v2y);
x_tmp=x(:,:,1)+alpha(:,:,1).*v2x; y_tmp=y(:,:,1)+alpha(:,:,1).*v2y;
beta(:,:,1)=(x_all2(:,:,1)-x_tmp).^2+(y_all2(:,:,1)-y_tmp).^2;

end

%at this point: beta contains the distance to the four sides, and this 
%is all we are going to need

for icur=1:4;
beta=circshift(beta,[0 0 1]); coeffs=circshift(coeffs,[0 0 1]);
coeffs(:,:,1)=beta(:,:,2).*beta(:,:,3);
end
coeffs=reshape(coeffs,[size(coeffs,1)*size(coeffs,2) size(coeffs,3)]);
coeffs=coeffs./(sum(coeffs,2)*ones(1,4));
coeffs=reshape(coeffs,[size(alpha,1) size(alpha,2) size(alpha,3)]);

warning on MATLAB:divideByZero;