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