function xy = uniform_in_hexagon01 ( n )

%*****************************************************************************80
%
%% uniform_in_hexagon01 samples uniformly from the regular unit hexagon.
%
%  Discussion:
%
%    The unit hexagon has center (0,0) and "radius" 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    19 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Output, real XY(2,N), the points.
%

%
%  Set basis vectors for each rhombus.
%
  v = [ ...
    -1.0, 0.5,            0.5,           -1.0; ...
     0.0, sqrt(3.0)/2.0, -sqrt(3.0)/2.0,  0.0 ];
%
%  Assign each point randomly to one of the rhombuses.
%
  rh = randi ( [ 1, 3 ], n, 1 );

  i = find ( rh == 1 );
  j = find ( rh == 2 );
  k = find ( rh == 3 );
%
%  Set barycentric coordinates of each point in its rhombus.
%
  xy = rand ( 2, n );
%
%  Convert to physical coordinates.
%
  xy(1:2,i) = v(1:2,1:2) * xy(1:2,i);
  xy(1:2,j) = v(1:2,2:3) * xy(1:2,j);
  xy(1:2,k) = v(1:2,3:4) * xy(1:2,k);

  return
end