gael/profilesMatlabProcessing/netcdf_ecco_GenericgridCoeffs_4points.m

%
%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;

1	gforget	1.1	%
2			%This function does a "near-bilinear interpolation".
3			%
4			%We have to work with distorted quadrilaterals rather than rectangles, but
5			%we use the distance to the quadrilateral sides as in bilinear interpolation:
6			%-> when the quadrilateral is a rectangle this actually is bilinear interpolation.
7			%-> for convex quadrilaterals it is close.
8			%-> for non-convex quadrilaterals !?
9			%
10			%Author: Gael Forget
11			%Date: June 14th 2007
12			%
13			function [coeffs]=netcdf_ecco_GenericgridCoeffs(x0,y0,x_all,y_all);
14
15			warning off MATLAB:divideByZero;
16
17			x=ones([length(x_all) size(x0)]); for icur=1:length(x_all); x(icur,:,:)=x0; end;
18			y=ones([length(x_all) size(x0)]); for icur=1:length(x_all); y(icur,:,:)=y0; end;
19
20			alpha=zeros(size(x)); beta=zeros(size(x)); gamma=zeros(size(x)); coeffs=zeros(size(x));
21
22			x_all2=zeros([length(x_all) size(x0,1) 1]); for icur=1:length(x_all); x_all2(icur,:,:)=x_all(icur); end;
23			y_all2=zeros([length(y_all) size(x0,1) 1]); for icur=1:length(y_all); y_all2(icur,:,:)=y_all(icur); end;
24
25			for icur=1:4;
26
27			x=circshift(x,[0 0 1]); y=circshift(y,[0 0 1]);
28			alpha=circshift(alpha,[0 0 1]); beta=circshift(beta,[0 0 1]);
29
30			v1x=x_all2(:,:)-x(:,:,1); v1y=y_all2(:,:,1)-y(:,:,1); v2x=x(:,:,2)-x(:,:,1); v2y=y(:,:,2)-y(:,:,1);
31			alpha(:,:,1)=(v1x.v2x+v1y.v2y)./(v2x.v2x+v2y.v2y);
32			x_tmp=x(:,:,1)+alpha(:,:,1).v2x; y_tmp=y(:,:,1)+alpha(:,:,1).v2y;
33			beta(:,:,1)=(x_all2(:,:,1)-x_tmp).^2+(y_all2(:,:,1)-y_tmp).^2;
34
35			end
36
37			%at this point: beta contains the distance to the four sides, and this
38			%is all we are going to need
39
40			for icur=1:4;
41			beta=circshift(beta,[0 0 1]); coeffs=circshift(coeffs,[0 0 1]);
42			coeffs(:,:,1)=beta(:,:,2).*beta(:,:,3);
43			end
44			coeffs=reshape(coeffs,[size(coeffs,1)*size(coeffs,2) size(coeffs,3)]);
45			coeffs=coeffs./(sum(coeffs,2)*ones(1,4));
46			coeffs=reshape(coeffs,[size(alpha,1) size(alpha,2) size(alpha,3)]);
47
48			warning on MATLAB:divideByZero;
49