function [zgnew,wg,rmserror,zr]=barnes(xx,yy,zg,x,y,z,zr,radius,ehalf,geomx,geomy)
% [zgnew,wg,rmserror,zr]=barnes(xx,yy,zg,x,y,z,zr,radius,ehalf,geomx,geomy)

zgnew=zeros(size(zg));
wg=zeros(size(zg));

% interpolate residual to grid
for i=1:length(xx)
  xge=xx(i);
  for j=1:length(yy)
    yge=yy(j);
    w=exp(-((x-xge)/radius*geomx).^2-((y-yge)/radius*geomy).^2);
    zgnew(j,i)=zg(j,i)+zr'*w/(sum(w)+ehalf);
    wg(j,i)=sum(w);
    end
  end

% interpolate from grid to data
tab=[[0;yy],[xx;zgnew]];
zest=zeros(size(z));
for k=1:length(z)
  zest(k)=table2(tab,y(k),x(k));
  end;
rmserror=sqrt(sum((z-zest).^2)/(length(z)-1));

% evaluate residual
zr=z-zest;
return