function [vv]=flowout(Z,i,j)
// the flow function from a point (i,j) in a domain Z
// output is an array of 3-tuples
// (i*,j*,f*) where f* is the fraction flowing in the neighbour (i*,j*)

xx=[];
ssum=0;

for ii=-1:2:1  // the X-neighbours
  if Z(i+ii,j)<Z(i,j);
    xx=[xx;(Z(i,j)-Z(i+ii,j)) i+ii j];
    ssum=ssum+(Z(i,j)-Z(i+ii,j));
   end;
end;

for jj=-1:2:1  // the Y-neighbours
  if Z(i,j+jj)<Z(i,j);
    xx=[xx;(Z(i,j)-Z(i,j+jj)) i j+jj];
    ssum=ssum+(Z(i,j)-Z(i,j+jj));
   end;
end;

for ii=-1:2:1  // the diagonal-neighbours
for jj=-1:2:1
  if Z(i+ii,j+jj)<Z(i,j);
    xx=[xx;(Z(i,j)-Z(i+ii,j+jj)) i+ii j+jj];
    ssum=ssum+(Z(i,j)-Z(i+ii,j+jj));
   end;
end;
end;
// ssum is the sum of all downward height differences

xx=gsort(xx,'lr');
// doesnt really matter, but just so that neighbours in decreasing order

vv=[xx(:,2) xx(:,3) xx(:,1)/ssum];

return;
endfunction;

