function h=lndfill(fem_grid_struct,color)
%LNDFILL draw land polygons for FEM domain
%   LNDFILL draws the land polygons generated by LNDMAKE
%   for a given FEM domain.  This is the "masking" routine
%   for the finite-difference based filled contouring package
%   FDCONT. The necessary fem_grid_struct fields for LNDFILL
%   are generated by LNDMAKE, and are contained in the
%   fem_grid_struct returned by LOADGRID, if the appropriate
%   files are located.
%
%   Input: fem_grid_struct - FEM domain grid structure (see LOADGRID)
%   	   color    - polygon face color
%
%   Output : h - handle to patch objects drawn
%
%   Call as: h=lndfill(fem_grid_struct,color);
%
%   Written by: Chris E. Naimie
%   Modified by: Brian O. Blanton for more general usage. (Jan 99)

if nargin<1 | nargin>2
   error('Incorrect number of input arguments to LNDFILL')
end
if nargout>1  
   error('Incorrect number of output arguments to LNDFILL')
end

% Verify incoming structure
if ~is_valid_struct(fem_grid_struct)
   error('    Argument to LNDFILL must be a valid fem_grid_struct.')
end

% Check to see if the .lnd and .lbe fields of fem_grid_struct
% exist and are not empty.
if ~isfield(fem_grid_struct,'lnd')
   error('    Land Node field not part of fem_grid_struct')
elseif ~isfield(fem_grid_struct,'lbe')
   error('    Land Element field not part of fem_grid_struct')
end
lnd=fem_grid_struct.lnd;
lbe=fem_grid_struct.lbe;
if isempty(lnd) 
   error('    Land Node field in fem_grid_struct is EMPTY')
elseif isempty(lbe) 
   error('    Land Element field in fem_grid_struct is EMPTY')
end

% If color is not passed in, set the color to
% draw the land polys to be the same as the gca 
% background color. 
if nargin == 1
   color=get(gca,'Color');
end

%
% Load data from .lnd and .lbe files
x=lnd(:,1);y=lnd(:,2);

%
% colorfill land and islands
i2=0;
count=0;
while i2 < length(lbe)
%
   i1=i2+1;
   i2=find(lbe(:,3)==lbe(i1,2));
   polygon=lbe(i1:i2,2);
   xp=x(polygon);yp=y(polygon);
   count=count+1;
   hf(count)=patch(xp,yp,color);
   set(hf(count),'EdgeColor','none')
end
if nargout==1
   h=hf;
end