function [bathynew,D] = smoobathy(elems,nodes,bathy)
% [bathynew, D] = smoobathy(meshele,meshnod,bathy);
% D is the smoothing matrix    
% bathynew(:,2) = D*bathy(:,2);

polys = elems;    
polys(:,1)=elems(:,4);   %   triangles close on themselves

[m,n]=size(polys);
c = ones(4*m,1);
c(:) = polys;
l = min(c):max(c);    % this will be harder. It leaves out unused nodes   
mc = min(c)-1;

xy=nodes(l,2:3);
i1 = polys(:,1)-mc;
j1 = polys(:,2)-mc;
i2 = polys(:,2)-mc;
j2 = polys(:,3)-mc;
i3 = polys(:,3)-mc;
j3 = polys(:,4)-mc;
i4 = polys(:,4)-mc;
j4 = polys(:,1)-mc;
S = sparse([i1;i2;i3],[j1;j2;j3],1);

%figure(1)
% spy(S)
%figure(2)
% gplot(S,xy)
%axis('equal')

%figure(3)
%Sr=symrcm(S);
%spy(S(Sr,Sr))  



%Create new bathy by clever smoothing
bathynew = bathy;

D = S;
[n,m]=size(S);

for rows = 1:n
l = find(S(rows,:)>0);
d = 1 ./sqrt((xy(rows,1)-xy(l,1)).^2 + (xy(rows,2)-xy(l,2)).^2);
d = d/sum(d);
l = [l ; rows];
d = [d; 1]/2;    % The present point  (rows,rows) gets half the weighting.
D(rows,l) = d';
end
bathynew(:,2) = D*bathy(:,2);



X=polys;
X(:)=xy(polys-mc,1);
X=X';
Y=polys;
Y(:)=xy(polys-mc,2);
Y=Y';
B = polys;
B(:) = bathy(polys-mc,2);
B=B';
figure(1)
clf
fill(X,Y,B)
title('before')
axis('equal')
B = polys;
B(:) = bathynew(polys-mc,2);
B=B';
figure(2)
clf
fill(X,Y,B)
title('after')
axis('equal')

return

% Reconstruct bathy from B and polys
% Go through all polys and set bathynew to at each node of poly
% If B and poly are complete then all of bathy is filled
for i = 1:length(polys(:,1)-mc)
k = polys(i,1)-mc;
bathynew(k,2) = Bnew(i,1);
k = polys(i,2)-mc;
bathynew(k,2) = Bnew(i,2);
k = polys(i,3)-mc;
bathynew(k,2) = Bnew(i,3);
k = polys(i,4)-mc;
bathynew(k,2) = Bnew(i,4);
end
disp('wrong way')

