
function [z2,ez2]=interp2_opt(x,y,z,x2,y2,corrdist,noise);

% function [z2,ez2]=interp2_opt(x,y,z,x2,y2,corrdist,noise);
% 
% Optimal interpolation to grid x2,y2, assuming that
% the observations z(x,y) are the sum of a true
% signal zt, with covariance 
%   var(z)*exp(-dist^2/corrdist^2)
% and white noise with covariance "noise".

meanz=mean(z);
z=z-meanz;

transflag=0;
if (length(x2)==size(x2,2))
    x2=x2';
    transflag=1;
end;

if ~exist('corrdist')
  corrdist=.5*mean(diff(x));
end;
if ~exist('noise')
  noise=0;
end;

x2_orig=x2; y2_orig=y2;

E=zeros(length(x)+length(x2), ...
    length(x)+length(x2));
E(1:length(x),1:length(x))=diag(ones(length(x),1));
b=[z;zeros(size(x2))];
if (length(noise)>1)
    noise=[noise;mean(noise)*ones(size(x2))];
end;
x2=[x;x2];y2=[y;y2];

% additional constraint: overall means equal
%E(size(E,1)+1,:)=1;
%b(length(b)+1)=length(x2)*sum(y)/length(y);

dX=abs(ones(length(x2),1)*x2' -x2*ones(1,length(x2)));
dY=abs(ones(length(y2),1)*y2' -y2*ones(1,length(y2)));
del2=dX.^2+dY.^2;
Rzz=var(z)*exp(-del2./(corrdist.^2));

if length(noise)==1
    Rnn=noise^2*diag(ones(length(x2),1));
else
    Rnn=diag(noise.^2);
end;
C=(Rzz*E'/(E*Rzz*E'+Rnn));
z2=C*b;

P=Rzz-C*E*Rzz;
ez2=sqrt(diag(P));

z2=z2(length(x)+1:length(z2));
ez2=ez2(length(x)+1:length(ez2));

if transflag
    z2=z2';
    ez2=ez2';
end;

z2=z2+meanz;

