function G = numgrid ( R, n )

%*****************************************************************************80
%
%% numgrid() numbers the grid points in a two dimensional region.
%
%  Discussion:
%
%    G = NUMGRID(R,N) numbers the points on an N-by-N grid in
%    the subregion of -1<=x<=1 and -1<=y<=1 determined by REGION.
%
%  Licensing:
%
%    Based on code that is copyright 1984-2015 The MathWorks, Inc. 
%
%  Modified:
%
%    30 August 2023
%
%  Author:
%
%    Original MATLAB version by the MathWorks.
%    This version by John Burkardt
%
%  Input:
%
%    character R: specifies the region:
%      'S' - the entire square.
%      'L' - the L-shaped domain made from 3/4 of the entire square.
%      'C' - like the 'L', but with a quarter circle in the 4-th square.
%      'D' - the unit disk.
%      'A' - an annulus.
%      'H' - a heart-shaped cardioid.
%      'B' - the exterior of a "Butterfly".
%
%    integer n: the number of rows and columns in the grid.
%
%  Output:
%
%    integer g(n,n): positions corresponding to unused nodes are 0.
%    Nodes included in the grid are numbered consecutively, starting
%    at the top left, proceeding downward, and then to the right.
%
  if ( n < 2 )
    error ( message ( 'numgrid(): Invalid number of points.' ) );
  end
%
%  Set the (x,y) coordinates of the nodes.
%
  x = ones ( n, 1 ) * linspace ( -1.0, +1.0, n );
  y = flipud ( x' );
%
%  Depending on the region, select those nodes which are active.
%
  if ( R == 'S' )
    G = (x > -1) & (x < 1) & (y > -1) & (y < 1);
  elseif ( R == 'L' )
    G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ( (x > 0) | (y > 0));
  elseif ( R == 'C' )
    G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ((x+1).^2 + (y+1).^2 > 1);
  elseif ( R == 'D' )
    G = x.^2 + y.^2 < 1;
  elseif ( R == 'A' )
    G = ( x.^2 + y.^2 < 1) & ( x.^2 + y.^2 > 1/3 );
  elseif ( R == 'H' )
    RHO = 0.75;
    SIGMA = 0.75;
    G = ( x.^2 + y.^2 ) .* ( x.^2 + y.^2 - SIGMA * y ) < RHO * x.^2;
  elseif ( R == 'B' )
    t = atan2 ( y, x );
    r = sqrt ( x.^2 + y.^2 );
    G = ( r >= sin ( 2.0 * t ) + 0.2 * sin ( 8.0 * t ) ) & ...
        (x > -1) & (x < 1) & (y > -1) & (y < 1);
  else
     error ( message ( 'numgrid(): Invalid region type.' ) );
  end
%
%  Convert from logical to integer.
%
  k = find ( G );
  G = zeros ( size ( G ) );
  G(k) = ( 1 : length ( k ) )';

  return
end

