function [coeffs]=netcdf_ecco_GenericgridCoeffs(x0,y0,x_all,y_all);

%This script is a two point interpolation scheme, which yields more continuous
%results than a one point "closest neighbor" interpolation scheme for example.
%
%Within a quadrilateral, the interpolated value is the projection of the 
%observed point on the closest side.
%
%As with the closest point scheme, there is no discontinuity through edges 
%separating two neighboring quadrilaterals. 
%

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]); gamma=circshift(gamma,[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;

gamma(:,:,1)=(v1x.*v1x+v1y.*v1y);

end

beta=reshape(beta,[size(beta,1)*size(beta,2) size(beta,3)]);
betamin=min(beta,[],2);
for icur=1:4;
beta=circshift(beta,[0 1]);
tmp1=find(beta(:,1)==betamin); 
beta(tmp1,1)=1; beta(tmp1,2:4)=0;
end
beta=reshape(beta,[size(alpha,1) size(alpha,2) size(alpha,3)]);


for icur=1:4;
alpha=circshift(alpha,[0 0 1]); beta=circshift(beta,[0 0 1]); coeffs=circshift(coeffs,[0 0 1]);
coeffs(:,:,1)=coeffs(:,:,1)+(1-alpha(:,:,1)).*beta(:,:,1);
coeffs(:,:,2)=coeffs(:,:,2)+alpha(:,:,1).*beta(:,:,1);
end

%revert to closest neighbor if very close, to avoid computer precision issue:
gamma=reshape(gamma,[size(gamma,1)*size(gamma,2) size(gamma,3)]); gamma=gamma./(max(gamma,[],2)*ones(1,4));
coeffs=reshape(coeffs,[size(coeffs,1)*size(coeffs,2) size(coeffs,3)]); 
for icur=1:4;
coeffs=circshift(coeffs,[0 0 1]); gamma=circshift(gamma,[0 0 1]);
tmp1=find(squeeze(gamma(:,1))>0.99);
coeffs(tmp1,1)=1; coeffs(tmp1,2:4)=0;
end
coeffs=reshape(coeffs,[size(alpha,1) size(alpha,2) size(alpha,3)]);