function s = polygon_sample ( nv, v, n )

%*****************************************************************************80
%
%% polygon_sample() uniformly samples a polygon.
%
%  Discussion:
%
%    The indexing of triangles was changed to (n-2,3), on 14 January 2025.
%    The indexing of polygon vertices was changed to (nv,2), 15 January 2025.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 January 2025
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NV, the number of vertices.
%
%    real V(NV,2), the polygon vertices in counterclockwise order.
%
%    integer N, the number of points to create.
%
%  Output:
%
%    real S(N,2), the points.
%

%
%  Triangulate the polygon.
%
  triangles = polygon_triangulate ( nv, v(1:nv,1), v(1:nv,2) );
%
%  Determine the areas of each triangle.
%
  for i = 1 : nv - 2
    area_triangle(i) = triangle_area ( ...
      v(triangles(i,1),1), v(triangles(i,1),2), ...
      v(triangles(i,2),1), v(triangles(i,2),2), ...
      v(triangles(i,3),1), v(triangles(i,3),2) );
  end
%
%  Normalize the areas.
%
  area_polygon = sum ( area_triangle(1:nv-2) );
  area_relative(1:nv-2) = area_triangle(1:nv-2) / area_polygon;
%
%  Replace each area by the sum of itself and all previous ones.
%
  area_cumulative(1) = area_relative(1);
  for i = 2 : nv - 2
    area_cumulative(i) = area_relative(i) + area_cumulative(i-1);
  end

  for j = 1 : n
%
%  Choose triangle I at random, based on areas.
%
    area_percent = rand ( 1, 1 );

    for k = 1 : nv - 2

      i = k;

      if ( area_percent <= area_cumulative(k) )
        break
      end

    end
%
%  Now choose a point at random in triangle I.
%
    r = rand ( 2, 1 );

    if ( 1.0 < sum ( r(1:2) ) )
      r(1:2) = 1.0 - r(1:2);
    end

    s(j,1:2) = ( 1.0 - r(1) - r(2) ) * v(triangles(i,1),1:2) ...
                     + r(1)          * v(triangles(i,2),1:2) ...
                            + r(2)   * v(triangles(i,3),1:2);
  end

  return
end

