function xy = sample_quad ( quad_xy, n )

%*****************************************************************************80
%
%% sample_quad() returns random points in a quadrilateral.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 February 2009
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real QUAD_XY(2,4), the coordinates of the nodes.
%
%    integer N, the number of points to sample.
%
%  Output:
%
%    real XY(2,N), the sample points.
%
  t1(1:2,1) = quad_xy(1:2,1);
  t1(1:2,2) = quad_xy(1:2,2);
  t1(1:2,3) = quad_xy(1:2,3);

  area1 = triangle_area ( t1 );

  t2(1:2,1) = quad_xy(1:2,3);
  t2(1:2,2) = quad_xy(1:2,4);
  t2(1:2,3) = quad_xy(1:2,1);

  area2 = triangle_area ( t2 );

  if ( area1 < 0.0 | area2 < 0.0 )

    t1(1:2,1) = quad_xy(1:2,2);
    t1(1:2,2) = quad_xy(1:2,3);
    t1(1:2,3) = quad_xy(1:2,4);

    area1 = triangle_area ( t1 );

    t2(1:2,1) = quad_xy(1:2,4);
    t2(1:2,2) = quad_xy(1:2,1);
    t2(1:2,3) = quad_xy(1:2,2);

    area2 = triangle_area ( t2 );

    if ( area1 < 0.0 | area2 < 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SAMPLE_QUAD - Fatal error!\n' );
      fprintf ( 1, '  The quadrilateral nodes seem to be listed in\n' );
      fprintf ( 1, '  the wrong order, or the quadrilateral is\n' );
      fprintf ( 1, '  degenerate.\n' );
      error ( 'SAMPLE_QUAD - Fatal error!' );
    end

  end

  area_total = area1 + area2;
%
%  Choose a triangle at random, weighted by the areas.
%  Then choose a point in that triangle.
%
  for i = 1 : n

    r = rand ( 1, 1 );

    if ( r * area_total < area1 )
      xy(1:2,i) = triangle_sample ( t1, 1 );
    else
      xy(1:2,i) = triangle_sample ( t2, 1 );
    end

  end

  return
end
