| 1 |
gforget |
1.1 |
function [coeffs]=netcdf_ecco_GenericgridCoeffs(x0,y0,x_all,y_all); |
| 2 |
|
|
|
| 3 |
|
|
%This script is a two point interpolation scheme, which yields more continuous |
| 4 |
|
|
%results than a one point "closest neighbor" interpolation scheme for example. |
| 5 |
|
|
% |
| 6 |
|
|
%Within a quadrilateral, the interpolated value is the projection of the |
| 7 |
|
|
%observed point on the closest side. |
| 8 |
|
|
% |
| 9 |
|
|
%As with the closest point scheme, there is no discontinuity through edges |
| 10 |
|
|
%separating two neighboring quadrilaterals. |
| 11 |
|
|
% |
| 12 |
|
|
|
| 13 |
|
|
x=ones([length(x_all) size(x0)]); for icur=1:length(x_all); x(icur,:,:)=x0; end; |
| 14 |
|
|
y=ones([length(x_all) size(x0)]); for icur=1:length(x_all); y(icur,:,:)=y0; end; |
| 15 |
|
|
|
| 16 |
|
|
alpha=zeros(size(x)); beta=zeros(size(x)); gamma=zeros(size(x)); coeffs=zeros(size(x)); |
| 17 |
|
|
|
| 18 |
|
|
x_all2=zeros([length(x_all) size(x0,1) 1]); for icur=1:length(x_all); x_all2(icur,:,:)=x_all(icur); end; |
| 19 |
|
|
y_all2=zeros([length(y_all) size(x0,1) 1]); for icur=1:length(y_all); y_all2(icur,:,:)=y_all(icur); end; |
| 20 |
|
|
|
| 21 |
|
|
for icur=1:4; |
| 22 |
|
|
|
| 23 |
|
|
x=circshift(x,[0 0 1]); y=circshift(y,[0 0 1]); |
| 24 |
|
|
alpha=circshift(alpha,[0 0 1]); beta=circshift(beta,[0 0 1]); gamma=circshift(gamma,[0 0 1]); |
| 25 |
|
|
|
| 26 |
|
|
v1x=x_all2(:,:)-x(:,:,1); v1y=y_all2(:,:,1)-y(:,:,1); v2x=x(:,:,2)-x(:,:,1); v2y=y(:,:,2)-y(:,:,1); |
| 27 |
|
|
alpha(:,:,1)=(v1x.*v2x+v1y.*v2y)./(v2x.*v2x+v2y.*v2y); |
| 28 |
|
|
x_tmp=x(:,:,1)+alpha(:,:,1).*v2x; y_tmp=y(:,:,1)+alpha(:,:,1).*v2y; |
| 29 |
|
|
beta(:,:,1)=(x_all2(:,:,1)-x_tmp).^2+(y_all2(:,:,1)-y_tmp).^2; |
| 30 |
|
|
|
| 31 |
|
|
gamma(:,:,1)=(v1x.*v1x+v1y.*v1y); |
| 32 |
|
|
|
| 33 |
|
|
end |
| 34 |
|
|
|
| 35 |
|
|
beta=reshape(beta,[size(beta,1)*size(beta,2) size(beta,3)]); |
| 36 |
|
|
betamin=min(beta,[],2); |
| 37 |
|
|
for icur=1:4; |
| 38 |
|
|
beta=circshift(beta,[0 1]); |
| 39 |
|
|
tmp1=find(beta(:,1)==betamin); |
| 40 |
|
|
beta(tmp1,1)=1; beta(tmp1,2:4)=0; |
| 41 |
|
|
end |
| 42 |
|
|
beta=reshape(beta,[size(alpha,1) size(alpha,2) size(alpha,3)]); |
| 43 |
|
|
|
| 44 |
|
|
|
| 45 |
|
|
for icur=1:4; |
| 46 |
|
|
alpha=circshift(alpha,[0 0 1]); beta=circshift(beta,[0 0 1]); coeffs=circshift(coeffs,[0 0 1]); |
| 47 |
|
|
coeffs(:,:,1)=coeffs(:,:,1)+(1-alpha(:,:,1)).*beta(:,:,1); |
| 48 |
|
|
coeffs(:,:,2)=coeffs(:,:,2)+alpha(:,:,1).*beta(:,:,1); |
| 49 |
|
|
end |
| 50 |
|
|
|
| 51 |
|
|
%revert to closest neighbor if very close, to avoid computer precision issue: |
| 52 |
|
|
gamma=reshape(gamma,[size(gamma,1)*size(gamma,2) size(gamma,3)]); gamma=gamma./(max(gamma,[],2)*ones(1,4)); |
| 53 |
|
|
coeffs=reshape(coeffs,[size(coeffs,1)*size(coeffs,2) size(coeffs,3)]); |
| 54 |
|
|
for icur=1:4; |
| 55 |
|
|
coeffs=circshift(coeffs,[0 0 1]); gamma=circshift(gamma,[0 0 1]); |
| 56 |
|
|
tmp1=find(squeeze(gamma(:,1))>0.99); |
| 57 |
|
|
coeffs(tmp1,1)=1; coeffs(tmp1,2:4)=0; |
| 58 |
|
|
end |
| 59 |
|
|
coeffs=reshape(coeffs,[size(alpha,1) size(alpha,2) size(alpha,3)]); |
| 60 |
|
|
|
| 61 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph