function geosurf( minlon, maxlon, minlat, maxlat, infile )
%
%    The function geosurf(minlon, maxlon, minlat, maxlat, infile)
%    makes a surface 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.
%    It redistributes the height data in a SQRT(N) x SQRT(N) 
%    matrix, so that they can be contoured.
%
%---------------------------------------------------------------- 
%    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;
%----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
%
% --- 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 ------
xx = [ xmin : xstep : xmax];
yy = [ ymin : ystep : ymax];
if (length(xx) ~= n), 
  ll = length(xx) - n;
  if ll == 2,
   x = xx(2:length(xx)-1);
  end
  if ll == 1,
   x = xx(1:length(xx)-1);
  end
 else
  x = xx;
end
if (length(yy) ~= m), 
  ll = length(yy) - m;
  if ll == 2,
   y = yy(2:length(yy)-1);
  end
  if ll == 1,
   y = yy(1:length(yy)-1);
  end
 else
  y = yy;
end
clear xx;
clear yy;
size(x)
size(y)
surf(x,y,Z)
