function riverdraw(efile, afile, vertangle, rotangle, threshold)
%
%  Written by: Glenn E. Moglen
%  
%  Script plots 3-D view of drainage basin with channel network shaded
%  using threshold area.  
%
%  efile:      file containing matrix of elevation values. (.dat extension assumed)
%  afile:      file containing matrix of cumulative area values.  (.dat extension assumed)
%  vertangle:  azimuth angle (0 => profile view, 90 => plan view)
%  rotangle:   angle relative to North-South (0 => looking north, 90 => looking west) 
%  threshold:  channel source area in pixels.
%
%  Note: efile and afile entries should be put in single quotes and 
%        should be entered without the assumed file extension of '.dat'.
%
%  Example: 
% 
%      >> riverdraw ('plumelev', 'plumarea', 70, 10, 100)
%
%  Matlab will plot 3-D image of elevations contained in 'plumelev.dat', and shade
%  the channel network at a threshold of 100 pixels based on the areas in
%  'plumarea.dat'.  The azimuth view angle is 70 degrees and rotational view
%  angle is 10 degrees.
%
   global elev; global area; global prop;
   eval(['load ', efile, '.dat'])
   eval(['load ', afile, '.dat'])
   e = eval(efile);
   a = eval(afile);
   [nx ny] = eval(['size(', efile, ')'])
   ctr = 0.0;
   minelev = 10000000;
   area_at_outlet = max(max(a))
   for i = 1:nx,
      for j = 1:ny,
         if (a(i,j) == 0),
            e(i,j) = 0;
         end 
      end 
   end 
   for i = 1:nx,
      for j = 1:ny,
         if (a(i, j) > 0 & e(i, j) < minelev),
            minelev = e(i, j);
         end
      end
   end 
   for i = 1:nx, 
      for j = 1:ny,
         if (e(i, j) > 0),
            e(i, j) = e(i, j) - minelev;
            if (a(i, j) > threshold),
               a(i, j) =  -1;
               ctr = ctr + 1;
            else
               a(i, j) = 1;
            end
         else
            a(i, j) = 0;
         end 
      end
   end
   colormap(gray);
   surf(e,a), view(rotangle, vertangle), caxis([-1.5 1]), grid
   maximum_elevation = max(max(e))
   Drainage_Density = ctr / area_at_outlet