function geomake5( minlon, maxlon, minlat, maxlat, infile )
%
%    The function geomake5(minlon, maxlon, minlat, maxlat, infile)
%    makes a contour plot of topographic data, between the
%    coordinates minlon < x < maxlon, minlat < y < maxlat,
%    reading from the file 'infile', which should contain 
%    (altitude - bathymetry) data, given in a form of a Nx3
%    matrix, where the first column is latitude, the second 
%    longitude, and the third height from the Sea Surface.
%
%    This subroutine is different from geomake in that IT ASSUMES
%    a grid size 5'x5' for the data (ETOPO5 data set)
%
%    The resulting contours are saved in a file called contfile.dat
%    containing the isobaths at -4000, -3000, -2000, -1000,
%    -500 and -200 meters, as well as in corresponding files
%    named dep###.dat, with data from a single isobath.
%
%    The best way to plot the resulting isobaths is use of 
%    the subroutine seaplot.
%
%    See also geosurf, geocont, map2, bwmap2.
%
%---------------------------------------------------------------- 
%    Written by Vassilis Zervakis, in January 1996, while at
%               AINCO - Interocean, Madrid, Spain.
%----------------------------------------------------------------
%
%----read the file----------------------
     filename = infile;
     nl = length(filename);
     eval(['load ' filename]);
     filer = eval('filename(1:nl-4)');
     X = eval(filer);
     clear infile ;
     clear filer ;
     llon = X(:,2);
     llat = X(:,1);
     dep = X(:,3);
     clear filename;
     clear X;
     length(llon)
%----decide how to create the matrix ----
% ---decide the limits of the final box
     if min(llat) < minlat,
        ymin = minlat;
      else
        ymin = min(llat);
     end
     if max(llat) > maxlat,
        ymax = maxlat;
      else
        ymax = max(llat);
     end
     if min(llon) < minlon,
        xmin = minlon;
      else
        xmin = min(llon);
     end
     if max(llon) > maxlon,
        xmax = maxlon;
      else
        xmax = max(llon);
     end
     xmin
     xmax
     ymin
     ymax
%
% --- Replacing the following lines:  KNOWING IT's a 5'x5' grid,
% --- define the stepsizes:
     xgrid = 5.;
     ygrid = 5.;
     xstep1 = xgrid / 60.;
     xstep2 = xgrid / 60.;
     ystep1 = xgrid / 60.;
     ystep2 = xgrid / 60.;
     xstep = xstep1;
     ystep = ystep1;
% --- do the same operation to find the size of the step --
% --- at the topmost latitude -----------------------------
%     y1 = ymin;
%    dymin = 0.1;
%    for m = 2 : length(llat) - 1,
%      dy = llat(m) - y1;
%      if (dy > 0.0001) & (dy < dymin) & (llat(m) < ymax),
%         dymin = dy;
%         mmin = m;
%      end
%    end
%    y2 = llat(mmin);
%%   ystep1 = y2 - y1
% --- do the same operation to find the size of the step --
% --- at the topmost latitude -----------------------------
%    yN = ymax;
%    dymin = 0.1;
%    for m = 1 : length(llat) - 1,
%      dy = yN - llat(m);
%      if (dy > 0.0001) & (dy < dymin) & (llat(m) > ymin),
%         dymin = dy;
%         mmax = m;
%      end
%    end
%    yN1 = llat(mmax);
%    ystep2 = yN - yN1
% --- do the same operation to find the size of the step --
% --- at the easternmost longitude-------------------------
%    x1 = xmin;
%    dxmin = 0.1;
%    for m = 1 : length(llon) - 1,
%      dx = llon(m) - x1;
%      if (dx > 0.0001) & (dx < dxmin) & (llon(m) < xmax),
%         dxmin = dx;
%         mmin = m;
%      end
%    end
%    x2 = llon(mmin);
%    xstep1 = x2 - x1
% --- do the same operation to find the size of the step --
% --- at the westernmost longitude-------------------------
%    xN = xmax;
%    dxmin = 0.1;
%    for m = 1 : length(llon),
%      dx = xN - llon(m);
%      if (dx>0.0001) & (dx<dxmin) & (llon(m) > xmin),
%         dxmin = dx;
%         mmax = m;
%      end
%    end
%    xN1 = llon(mmax);
%    llon(mmax)
%    xmin
%    xstep2 = xN - xN1
%
% ---Compute a probable size for the matrix: -----------------
     xav = (xmax - xmin) / (0.5 * (xstep1 + xstep2));
     yav = (ymax - ymin) / (0.5 * (ystep1 + ystep2));
     m = round(yav);
     n = round(xav);
%
% ---start (at last) creating the desired matrix of eights ---
% ------------------------------------------------------------
%
     for k = 1:length(llat),
       if (xmin<llon(k))&(xmax>llon(k))&(ymin<llat(k))&(ymax>llat(k)),
%  --select column from longitude -----------------------------------
        if xstep2 ~= xstep1,
          xstep=((xstep2-xstep1)/(xmax-xmin))*(llon(k)-xmin)+xstep1;
         else
          xstep = xstep2;
        end
        kj = round((llon(k) - xmin) / xstep);
%  --select row from latitude----------------------------------------
        if ystep2 ~= ystep1,
          ystep=((ystep2-ystep1)/(ymax-ymin))*(llat(k)-ymin)+ystep1;
         else
          ystep = ystep2;
        end
        ki = m - round((ymax - llat(k)) / ystep);
%  --fill the matrix at last -------------------------------
        Z(ki,kj) = dep(k);
       end
     end
%
[m,n]=size(Z);
clear dep;
clear llon;
clear llat;
%  ---make the x vector to help plotting ------
for k = 1:m, 
  y(k) = ymin + (k-1)*ystep;
end
for k = 1:n,
  x(k) = xmin + (k-1)*xstep;
end
D = contour(x,y,Z,[-4000 -3000 -2000 -1000 -500 -200]);
clear Z;
clear x;
clear y;
E = D';
save contfile.dat E -ascii
% -----save the -4000 contour line -------------
pairs1 = D(2,1);
E = D(:,1:pairs1+1)';
dd = -D(1,1);
depth = num2str(dd);
newfile = ['dep' depth '.dat'];
pairs1 = pairs1 + 1;
while dd == -D(1, pairs1+1),
   pairs2 = D(2,pairs1 + 1);
   E = [ E ; D(:,pairs1 + 1:pairs2+1)'];
   pairs1 = pairs1 + 1;
end
eval(['save ' newfile ' E -ascii'])
clear E ;
% -----save the -3000 contour line -------------
pairs2 = D(2,pairs1+1) ;
E = D(:,pairs1+2:pairs2+pairs1+1)';
dd = -D(1,pairs1+1) ;
depth = num2str(dd);
newfile = ['dep' depth '.dat'];
pairs1 = pairs1 + pairs2 + 1;
while dd == -D(1, pairs1+1),
   pairs2 = D(2,pairs1 + 1);
   mm = pairs1 + 1;
   E = [ E ; D(:,mm:mm+pairs2)'];
   pairs1 = pairs1 + pairs2 + 1;
end
eval(['save ' newfile  ' E -ascii']) ;
clear E ;
% -----save the -2000 contour line -------------
pairs3 = D(2,pairs1+1);
E = D(:,pairs1+2:pairs3+1)';
dd = -D(1,pairs1+1);
depth = num2str(dd);
newfile = ['dep' depth '.dat'];
pairs1 = pairs1 + pairs3 + 1;
while dd == -D(1, pairs1+1),
   pairs3 = D(2,pairs1 + 1);
   mm = pairs1 + 1;
   E = [ E ; D(:,mm:mm+pairs3)'];
   pairs1 = pairs1 + pairs3 + 1;
end
eval(['save ' newfile  ' E -ascii']) ;
clear E ;
% -----save the -1000 contour line -------------
pairs4 = D(2,pairs1+1);
E = D(:,pairs1+2:pairs4+1)';
dd = -D(1,pairs1+1);
depth = num2str(dd);
newfile = ['dep' depth '.dat'];
pairs1 = pairs1 + pairs4 + 1;
while dd == -D(1, pairs1+1),
   pairs4 = D(2,pairs1 + 1);
   mm = pairs1 + 1;
   E = [ E ; D(:,mm:mm+pairs4)'];
   pairs1 = pairs1 + pairs4 + 1;
end
eval(['save ' newfile  ' E -ascii']) ;
clear E ;
% -----save the -500 contour line -------------
pairs5 = D(2,pairs1+1);
E = D(:,pairs1+2:pairs5+1)';
dd = -D(1,pairs1+1);
depth = num2str(dd);
newfile = ['dep' depth '.dat'] ;
pairs1 = pairs1 + pairs5 + 1;
while dd == -D(1, pairs1+1),
   pairs5 = D(2,pairs1 + 1);
   mm = pairs1 + 1;
   E = [ E ; D(:,mm:mm+pairs5)'];
   pairs1 = pairs1 + pairs5 + 1;
end
eval(['save ' newfile  ' E -ascii']) ;
clear E ;
% -----save the -200 contour line -------------
pairs6 = D(2,pairs1+1) ;
E = D(:,pairs1+2:pairs6+1)';
dd = -D(1,pairs1+1) ;
depth = num2str(dd);
newfile = ['dep' depth '.dat'] ;
pairs1 = pairs1 + pairs6 + 1 ;
while (dd == -D(1, pairs1+1))&(pairs1 < length(D(2,:) - 1)),
   pairs6 = D(2,pairs1 + 1);
   mm = pairs1 + 1;
   E = [ E ; D(:,mm:mm+pairs6)'];
   pairs1 = pairs1 + pairs6 + 1;
end
eval(['save ' newfile  ' E -ascii']) ;
clear E ;
% -----save the -0 contour line -------------
% pairs7 = D(2,pairs1+1)
% E = D(:,pairs1+2:pairs7+1)';
% dd = -D(1,pairs1+1);
% depth = num2str(dd)
% newfile = ['dep' depth '.dat']
% pairs1 = pairs1 + pairs7 + 1
% while dd == -D(1, pairs1+1),
%    pairs7 = D(2,pairs1 + 1);
%    mm = pairs1 + 1;
%    E = [ E ; D(:,mm:mm+pairs7)'];
%    pairs1 = pairs1 + pairs7 + 1
% end
% eval(['save ' newfile  ' E -ascii'])
% clear E
