function geofill( minlon, maxlon, minlat, maxlat, coastfile, bathyfile )
%
%    The function geofill(minlon, maxlon, minlat, maxlat, 
%                         coastfile,bathyfile)
%    makes a plot of black (filled-in) land and isobaths between the
%    coordinates minlon < x < maxlon, minlat < y < maxlat. The coast
%    line data should be read from the file 'coastfile', while the
%    bathymetric data (isobaths) from the file 'bathyfile', which 
%    if it doesn't exist, can be created by running my MATLAB sub-
%    routine 'geomake'.  This subroutine is MUCH faster than geosurf,
%    or geocont, and (on the contrary to the forementioned) can be
%    superimposed by float tracks.
%
%    See also geocont, geosurf, geomake
%
%---------------------------------------------------------------- 
%    Written by Vassilis Zervakis, in January 1996, while at
%               AINCO - Interocean, Madrid, Spain.
%----------------------------------------------------------------
%
%----read and plot the coastline file ----------
%    figure
     coastfile;
     eval('map2(minlon, maxlon, minlat, maxlat, coastfile)');
     hold on;
%----read the bathymetry file ----------
%filepath = 'c:\geomaps\bathy\';
filenm = num2str(bathyfile);
n = findstr('.',filenm);
dirs = findstr('\',filenm);
nn = length(dirs);
m = dirs(nn);
filebody = filenm(m+1:n-1);
eval(['load ',filenm]);
data = eval(filebody);
depth = data(1,1);
points= data(1,2);
x = data(2:points+1,1);
y = data(2:points+1,2);
plot(x,y,'b');
hold on;
clear x;
clear y;
points_read = points + 1;
fpoints = 0.;
xf = [];
yf = [];
x = [];
POINT1 = 'A';
while points_read < length(data(:,1)),
 m = points_read+1;
 depth = data(m,1);
 points = data(m,2);
 points_read = points_read + points + 1;
 for kk = m+1:points_read,
   x1 = data(kk,1);
   y1 = data(kk,2);
   if (x1>=minlon)&(x1<=maxlon)&(y1>=minlat)&(y1<=maxlat),
      x(kk) = x1;
      y(kk) = y1;
   end
 end 
 if (depth == -4000)&(x ~= []),
  nn = length(x);
  x(1)
  x(nn)
  y(1)
  y(nn) 
  if (abs(x(1)-x(nn))+abs(y(1)-y(nn))) < 0.05
     fill(x,y,[0 0.25 0.25]);
   else
      POINT1
      if (abs(x(1) - minlon)) < 0.05,
         xf = [ xf; minlon; x];
         yf = [ yf; y(1)  ; y];
      end
      if (abs(x(1) - maxlon)) < 0.05,
         xf = [ xf; maxlon; x];
         yf = [ yf; y(1)  ; y];
      end
      if (abs(x(nn) - minlon)) < 0.05,
         xf = [ xf; x; minlon];
         yf = [ yf; y; y(nn) ];
      end
      if (abs(x(nn) - maxlon)) < 0.05,
         xf = [ xf; x; maxlon];
         yf = [ yf; y;  y(nn)];
      end
      if (abs(y(1) - minlat)) < 0.05,
         yf = [ yf; minlat; y];
         xf = [ xf; x(1)  ; x];
      end
      if (abs(y(1) - maxlat)) < 0.05,
         yf = [ yf; maxlat; y];
         xf = [ xf; x(1)  ; x];
      end
      if (abs(y(nn) - minlat)) < 0.05,
         yf = [ yf; y; minlat];
         xf = [ xf; x; x(nn)];
      end
      if (abs(y(nn) - maxlat)) < 0.05,
         yf = [ yf; y; maxlat];
         xf = [ xf; x; x(nn) ];
      end
      if fpoints == 0,
        fpoints = points + 1;
       else
        fpoints = fpoints + points;
      end
   end
 end
 hold on;
 clear x;
 clear y;
end
fill(xf,yf,[0 0.55 0.55]) ;
   
