function s = polygon_sample ( nv, v, n )

%*****************************************************************************80
%
%% polygon_sample() uniformly samples a polygon.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NV, the number of vertices.
%
%    real V(2,NV), the vertices of the polygon, listed in
%    counterclockwise order.
%
%    integer N, the number of points to create.
%
%  Output:
%
%    real S(N,2), the points.
%
  s = zeros ( n, 2 );
%
%  Triangulate the polygon.
%
  triangles = polygon_triangulate ( nv, v(1,1:nv), v(2,1:nv) );
%
%  Determine the areas of each triangle.
%
  for i = 1 : nv - 2
    area_triangle(i) = triangle_area ( ...
      v(1,triangles(1,i)), v(2,triangles(1,i)), ...
      v(1,triangles(2,i)), v(2,triangles(2,i)), ...
      v(1,triangles(3,i)), v(2,triangles(3,i)) );
  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 ( );

    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(1:2,triangles(1,i))' ...
                     + r(1)          * v(1:2,triangles(2,i))' ...
                            + r(2)   * v(1:2,triangles(3,i))';
  end

  return
end

