function [ZZ,locn]=drain(Z,loc,i,j,f)
// finds the net flowout or local minimas where flow from the point (i,j) 
// will go
//
// ZZ stores the above flowouts and minimas
// locn stores the paths

[m n]=size(Z);
ZZ=zeros(m,n);

if i==1|i==m|j==1|j==n  // if on the border of domain then over
  ZZ(i,j)=f;
  locn=loc;
  return;
end;

vv=flowout(Z,i,j); // compute the flowouts from (i,j)

[k,ig]=size(vv);
if k==0     // local minima
  ZZ(i,j)=f;
  locn=loc;
  return;
end;

if f<0.01    // if flow is small, then flow moves without bifurcating
  if k>0
    loc(vv(1,1),vv(1,2))=f; // uses that largest is first
    [ZZ,locn]=drain(Z,loc,vv(1,1),vv(1,2),f);
  end;
  return;
end;

for ii=1:1:k  // for all the neighbours 
  linn=zeros(m,n);
  linn(vv(ii,1),vv(ii,2))=vv(ii,3)*f; // distribute the flows
  [Zx,lx]=drain(Z,linn,vv(ii,1),vv(ii,2),vv(ii,3)*f); // recurse
  ZZ=ZZ+Zx;    // do the tally
  loc=loc+lx;
end;

locn=loc;
return;


endfunction;
