function program_04 ( )

%*****************************************************************************80
%
%% program_04() demonstrates sampling of a triangle.
%
%  Discussion:
%
%    The program
%    * reads a triangle T (defined by three points),
%    * reads "1" for bad scheme, "2" for good scheme.
%    * reads N, the number of random values to generate;
%    * it then computes N random points in the triangle;
%    * it determines N1, N2, and N3, the number of points
%      that ended up in subtriangles mbc, amc, and abm,
%      where "m" is the centroid.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 February 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'program_04(): Uniform Sampling\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Define a triangle T:\n' );

  t_v1 = [   0.0, 0.0 ]; 
  t_v2 = [ 100.0, 0.0 ]; 
  t_v3 = [   1.0, 1.0 ]; 
%
%  Get the scheme to use.
%
  for scheme = 1 : 2

    if ( scheme == 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Use bad sampling scheme.\n' );
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Use good sampling scheme.\n' );
    end
%
%  Get the number of values to generate.
%
    n = 1000;

    n1 = 0;
    n2 = 0;
    n3 = 0;

    for i = 1 : n

      if ( scheme == 1 )
        r = rand();
        s = rand();
        t = rand();
        xi1 = r / ( r + s + t );
        xi2 = s / ( r + s + t );
        xi3 = t / ( r + s + t );
      elseif ( scheme == 2 )
        r = rand();
        s = rand();
        xi1 =   1.0       - sqrt ( s );
        xi2 = ( 1.0 - r ) * sqrt ( s );
        xi3 =         r   * sqrt ( s );
      end

      if ( xi1 < xi2 & xi1 < xi3 )
        n1 = n1 + 1;
      elseif ( xi2 < xi1 & xi2 < xi3 )
        n2 = n2 + 1;
      elseif ( xi3 < xi1 & xi3 < xi2 )
        n3 = n3 + 1;
      end

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N1 = %8d  %% = %8.4f\n', n1, 100 * n1 / n );
    fprintf ( 1, '  N2 = %8d  %% = %8.4f\n', n2, 100 * n2 / n );
    fprintf ( 1, '  N3 = %8d  %% = %8.4f\n', n3, 100 * n3 / n );
    fprintf ( 1, '  N  = %8d  %% = %8.4f\n', n, 100 );
  
  end

  return
end
