function chebyshev_test ( m1, m2, n )

%*****************************************************************************80
%
%% chebyshev_test() looks at [-1,+1]x[-1,1] with Chebyshev PDF.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M1, M2, the number of rectangles in X and Y directions.
%
%    integer N, the number of sample points to be generated.
%
  x_min = -1.0;
  x_max = +1.0;
  y_min = -1.0;
  y_max = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'chebyshev_pdf_test():\n' );
  fprintf ( 1, '  Region [-1,+1]x[0,1].\n' );
  fprintf ( 1, '  PDF(X,Y)=1/sqrt(1-X^2)/sqrt(1-Y^2)\n' );
  fprintf ( 1, '  Use M1xM2=%dx%d grid of rectangles\n', m1, m2 )
  fprintf ( 1, '  Generate %d samples\n', n );
%
%  XB and YB are the 1D break points.
%
  xb = r8vec_cheby_extreme  ( m1 + 1, x_min, x_max );
  yb = r8vec_cheby_extreme  ( m2 + 1, y_min, y_max );
%
%  XC and YC are the 1D center points.
%
  xc = r8vec_cheby_zero ( m1, x_min, x_max );
  yc = r8vec_cheby_zero ( m2, y_min, y_max );
%
%  XC_MAT and YC_MAT are the M1 * M2 coordinates of rectangle centers.
%
  xc_mat = zeros ( m1, m2 );
  for i = 1 : m1
    xc_mat(i,1:m2) = xc(i);
  end

  yc_mat = zeros ( m1, m2 );
  for j = 1 : m2;
    yc_mat(1:m1,j) = yc(j);
  end
%
%  Evaluate the PDF at each rectangle center.
%
  pdf_mat(1:m1,1:m2) = 1.0 ./ sqrt ( 1.0 - xc_mat(1:m1,1:m2) .^ 2 ) ...
                           ./ sqrt ( 1.0 - yc_mat(1:m1,1:m2) .^ 2 );
%
%  Compute the discrete integral of PDF_MAT.
%
  integral = 0.0;
  for j = 1 : m2
    for i = 1 : m1
      rectangle = ( xb(i+1) - xb(i) ) * ( yb(j+1) - yb(j) );
      integral = integral + pdf_mat(i,j) * rectangle;
      cdf_mat(i,j) = integral;
    end
  end
%
%  Normalized the CDF.
%
  cdf_mat(1:m1,1:m2) = cdf_mat(1:m1,1:m2) / integral;
%
%  Select N values U uniformly at random in [0,1].
%
  u = rand ( n, 1 );
%
%  The CDF maps each U to a rectangle, within which we can
%  select a point (X,Y) uniformly at random.
%
  s = discrete_cdf_to_xy ( m1, m2, cdf_mat, xb, yb, n, u );
%
%  Write data to a file for examination, plotting, or analysis.
%
  filename = 'chebyshev.txt';
  r8mat_write ( filename, 2, n, s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote sample data to file "%s".\n', filename );
%
%  Plot data.
%
  figure ( );
  filename = 'chebyshev.png';
  plot ( s(1,:), s(2,:), '.', 'Markersize', 10 );
  axis ( 'equal' );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( 'PDF=1/(1-x^2)/(1-y^2)' );
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );

  return
end
