function cg = annulus_grid ( n, r1, r2, c, ng )

%*****************************************************************************80
%
%% annulus_grid() computes grid points inside an annulus.
%
%  Discussion:
%
%    The grid is defined by specifying the center of the annulus,
%    the inner and outer radius, and the number of subintervals N 
%    into which the horizontal radial segment should be divided.  
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 November 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of subintervals.
%
%    real R1, R2, the inner and outer radius of the annulus.
%
%    real C(2), the coordinates of the center of the annulus.
%
%    integer NG, the number of grid points, as determined by annulus_grid_count.
%
%  Output:
%
%    real CG(NG,2), the grid points inside the disk.
%
  cg = zeros ( ng, 2 );

  p = 0;

  for j = 0 : n
    i = 0;
    x = c(1);
    y = c(2) + r2 * 2 * j / ( 2 * n + 1 );
    if ( r1 * r1 <= ( x - c(1) ).^2 + ( y - c(2) ).^2 )
      p = p + 1;
      cg(p,1:2) = [ x,            y ];
      if ( 0.0 < r1 || 0 < j )
        p = p + 1;
        cg(p,1:2) = [ x, 2 * c(2) - y ];
      end
    end

    while ( true )

      i = i + 1;
      x = c(1) + r2 * 2 * i / ( 2 * n + 1 );

      if ( r1 * r1 <= ( x - c(1) ).^2 + ( y - c(2) ).^2 )
        if ( r2 * r2 < ( x - c(1) ).^2 + ( y - c(2) ).^2 )
          break
        end
        p = p + 1;
        cg(p,1:2) = [            x,            y ];
        p = p + 1;
        cg(p,1:2) = [ 2 * c(1) - x ,           y ];
        if ( 0 < j )
          p = p + 1;
          cg(p,1:2) = [            x, 2 * c(2) - y ];
          p = p + 1;
          cg(p,1:2) = [ 2 * c(1) - x , 2 * c(2) - y ];
        end
      end
    end
  end

  return
end
