function [xy_grid,xgrid,ygrid] = rect_grid ( datfile, coastfile, var, extra_coasts )
% RECT_GRID:  Constructs rectangular grid from data.
%
% Usage:
%      [data_grid,xgrid,ygrid] = rect_grid ( datfile, coastfile, var )

% John Evans (jevans@sakhalin.er.usgs.gov)
%

if ( nargin < 3 )
  help make_grid;
end


eval ( sprintf ( 'load %s', datfile ) );
eval ( sprintf ( 'load %s', coastfile ) );

% plot the coastline
if (~exist('coast'))
  fprintf (1, 'No coastline variable in coastfile.\n' );
  return;
else
%  plot(coast(:,1),coast(:,2), '-r');
end


% determine min and max extents of the grid from the x,y data
xmin = floor(min(x))*0.95;
xmax = ceil(max(x))*1.05;
ymin = floor(min(y))*0.95;
ymax = ceil(max(y))*1.05;


% determine grid size.  Assume that the spacing is 3 klicks.
xsize = ceil((xmax-xmin)/3.0);
ysize = ceil((ymax-ymin)/3.0);
data_grid = zeros(ysize,xsize);
grid_bndry = [xmin ymin; xmax ymin; xmin ymax; xmax ymax];



% put up the grid
domain_grid_fig = figure ( 'name', 'rect_grid:  grid domain' );

xax = [xmin:3.0:xmin+3.0*xsize];
yax = [ymin:3.0:ymin+3.0*ysize];
[XAX,YAX] = meshgrid(xax,yax);
zax = 0*XAX;


% for each separate coastline, set all grid points contained within equal to
% -Inf
coastline_breaks = find(isnan(coast(:,1)));

% make sure we get the last coastline
n = length(coastline_breaks);
coastline_breaks(n+1) = length(coast(:,1)) + 1;

endi = 0;
%for i = 1:length(coastline_breaks)
for i = 1:1

  % extract the ith coastline
  starti = endi+1;
  endi = coastline_breaks(i);
  coasti = coast(starti:10:endi-1,:);
  
  n = length(coasti(:,1));
  coasti(n+1,:) = coasti(1,:);

  
  % check to see if the coastline intersects the grid boundary.  If not,
  % don't bother checking for whether or not any of the grid cells are
  % contained by the coastline.  This is sort of the inverse problem, but
  % this way if there is a small island just outside of the grid, we don't
  % have to check every friggin' grid cell to see if it's inside the island.
  inside_out_coast = insider ( coasti(:,1), coasti(:,2), ...
                               grid_bndry(:,1), grid_bndry(:,2) );
  if ( find(inside_out_coast) )

    % need to do each column of the grid.  Can't do it all at once because of
    % memory problems.
    for j = 1:length(xax)
  
      inside_or_out_indices = insider ( XAX(:,j), YAX(:,j), ...
	                                coasti(:,1), coasti(:,2) );
      inside = find ( inside_or_out_indices );
      data_grid(inside,j) = -Inf * ones(size(inside'));
      zax(inside,j) =  -Inf * ones(size(inside'));
  
    end %for j

  end %if
  

  
end  % for i



% Set all data cells that are not -Inf to Inf.  These are the cells to grid.
ifind = find(data_grid ~= -Inf);
data_grid(ifind) = Inf * ones(size(ifind));



mesh(XAX,YAX,zax);
hold on;
set(gca,'view',[0 90]);
plot ( x, y,  '*g' );
plot(coast(:,1),coast(:,2), '-c');


% Now fill in the data that we have.
for i = 1:length(x)
  xind = 1 + floor((x(i) - xmin)/3.0);
  yind = 1 + floor((y(i) - ymin)/3.0);
  data_grid(yind,xind) = eval ( sprintf('%s(%i)', var, i) );
end


xy_grid = data_grid;
xgrid = XAX;
ygrid = YAX;

return;
