function lndmake(fem_grid_struct,belfilename)
%LNDMAKE create land masking information for FD CONTOURING
%   LNDMAKE creates a polygon list from existing .bel  and the node
%   lists from a fem_grid_struct. This polygon list can be used with 
%   the FD-based contouring and vector plotting functions (FDCONT) 
%   to mask out the land.
%
%   The .bel file must have been generated by the OPNML boundary
%   code generator GENBEL, since this code orders the boundary 
%   elements in order, with islands connected and last in the file.
%   Other codes (CONVCODES, fortran, Dartmouth) also output .bel
%   information in the needed order.  Work is ongoing to remove
%   this requirement.  If a .bel file needs ordering, load the
%   .bel file into GENBEL and directly  output it to a different
%   filename.
%
%   LNDMAKE outputs new files for the land nodes and segments.
%   The new files are <gridname>.lnd and <gridname>.lbe, and are
%   output to the current working directory.  Thus, LNDMAKE need
%   only be run once per domain.
%
%   Input: fem_grid_struct - FEM domain structure
%   	   belfilename - FEM domain .bel file
%
%   Output: NONE (See above)
%
%   Call as: >> lndmake(fem_grid_struct,belfilename)
%
%   Written by: Chris E. Naimie
%   Modified by: Brian O. Blanton to more general usage. (Jan 99)

%-----------------------------------------------------------------------

if nargin ~=2
   error('    Incorrect number of input arguments to LNDMAKE');
end

if ~is_valid_struct(fem_grid_struct)
   error('    Argument to LNDMAKE must be a valid fem_grid_struct.')
end

x=fem_grid_struct.x;
y=fem_grid_struct.y;

% Read boundary elements
%
[bel,gridname]=read_bel(belfilename);
%
% Outer boundary
%
% find last non-island boundary segment
notsea=find(bel(:,5)~=2);
lastnotsea=max(notsea);
sea=find(bel(:,5)==2); firstsea=min(sea);
if(firstsea<lastnotsea)
  sea=lastnotsea+1:length(bel);
end
% find last non-island boundary segment
sea=find(bel(:,5)==2); 
% do a recurrency test within the 2's
[elts,nelt]=count(bel(sea,2:3)); 
% islands within outer boundaries have last boundary elements with 
% only one occurency within the bel file
oones=find(nelt==1); oones=reshape(oones,2,length(oones)/2)';
beg=sea(oones(:,1)); en=oones(:,2);

sea(elts(beg):elts(en))=[];
%
if ~isempty(sea)
   polygon=bel(1:sea(1)-1,2);
else
   polygon=bel(:,2);        % no islands
end

xp=x(polygon);yp=y(polygon);
%
xrange=max(xp)-min(xp);yrange=max(yp)-min(yp);
ni=length(xp)+1; xp(ni)=xp(1);              yp(ni)=yp(1);
ni=length(xp)+1; xp(ni)=min(xp)-0.1*xrange; yp(ni)=yp(ni-1);
ni=length(xp)+1; xp(ni)=xp(ni-1);           yp(ni)=max(yp)+0.1*yrange;
ni=length(xp)+1; xp(ni)=max(xp)+0.1*xrange; yp(ni)=yp(ni-1);
ni=length(xp)+1; xp(ni)=xp(ni-1);           yp(ni)=min(yp)-0.1*yrange;
ni=length(xp)+1; xp(ni)=xp(1)  -0.1*xrange; yp(ni)=yp(ni-1);
ni=length(xp)+1; xp(ni)=xp(1)  -0.1*xrange; yp(ni)=yp(1);
%
delete(findobj(0,'Tag','LNDMAKE fig'))
figure('MenuBar','none','Name','LNDMAKE Output','Tag','LNDMAKE fig')
title(['LNDMAKE on ' gridname])
whitebg('w');
hp=line(xp,yp,'Color','b','LineStyle','-','LineWidth',2);
axis('tight')
drawnow

%
gridname=fem_grid_struct.name;
lndname=[gridname,'.lnd'];
lbename=[gridname,'.lbe'];
fnod=fopen(lndname,'w');
fbel=fopen(lbename,'w');
ione=1;
for i=1:length(xp)-1
   fprintf(fnod,'%i %f %f\n',i,xp(i),yp(i));
   fprintf(fbel,'%i %i %i\n',i,i,i+1);
end
i=length(xp);
fprintf(fnod,'%i %f %f\n',i,xp(i),yp(i));
fprintf(fbel,'%i %i %i\n',i,i,ione);

if ~isempty(sea)
   %
   % first island
   %
   clear polygon;
   i1=sea(1);
   i2=find(bel(:,3)==bel(i1,2));
   polygon=bel(i1:i2,2);
   xp=x(polygon);yp=y(polygon);
   hp=line(xp,yp,'Color','r','LineStyle','-');
   drawnow
   %
   ione=i+1;
   for ii=1:length(xp)-1
      i=i+1;
      fprintf(fnod,'%i %f %f\n',i,xp(ii),yp(ii));
      fprintf(fbel,'%i %i %i\n',i,i,i+1);
   end
   i=i+1;
   fprintf(fnod,'%i %f %f\n',i,xp(length(xp)),yp(length(xp)));
   fprintf(fbel,'%i %i %i\n',i,i,ione);
   %
   % remaining islands
   %
   while i2 < length(bel)
   %
      i1=i2+1;
      i2=find(bel(:,3)==bel(i1,2));
      polygon=bel(i1:i2,2);
      xp=x(polygon);yp=y(polygon);
      hp=line(xp,yp,'Color','r','LineStyle','-');
      drawnow
   %
      ione=i+1;
      for ii=1:length(xp)-1
   	 i=i+1;
   	 fprintf(fnod,'%i %f %f\n',i,xp(ii),yp(ii));
   	 fprintf(fbel,'%i %i %i\n',i,i,i+1);
      end
      i=i+1;
      fprintf(fnod,'%i %f %f\n',i,xp(length(xp)),yp(length(xp)));
      fprintf(fbel,'%i %i %i\n',i,i,ione);
   end
end
fclose(fnod);
fclose(fbel);

disp([lndname ' output to ' pwd])
disp([lbename ' output to ' pwd])