function s = cvt_1d_sample ( n )

%*****************************************************************************80
%
%% cvt_1d_sample() samples a PDF for a 1D CVT problem.
%
%  Discussion:
%
%    The Chebyshev points X have an asymptotic density
%
%      rho(x) = 1/pi * 1/sqrt(1-x^2)
%
%    Denote by H(I) the size of any of the disjoint intervals
%    containing X(I), and note that it there is a sort of 
%    equidistribution property so that
%
%      h(x(i)) proportional to 1 / rho(x(i))
%
%    It is known that if points X2 forming a CVT in 1D are generated 
%    according to a density rho2(x2), and H2(I) represents the diameter of the 
%    CVT interval, that we expect
%
%      h2(x2(i)) proportional to 1 / rho2(x2(i))^3
%
%    Therefore, we expect that if we can use the density
%
%      rho2(x) = rho(x)^(1/3)
%
%    that the CVT points will tend to the locations of the Chebyshev points.
%
%
%    As it turns out, the computation of the CVT requires sampling with
%    respect to the density rho, which in turn requires computing the
%    incomplete integral or CDF(x) of the density, and then inverting to
%    get X(CDF).  For the given density function above, Mathematica
%    is unable to determine an inverse function.
%
%    Therefore, to compute X(CDF), we must try using rejection sampling.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of points to compute.
%
%  Output:
%
%    real S(N), the sample points.
%
  s = zeros ( n, 1 );
%
%  OK, first we do non-vectorized.
%
  i = 0;

  while ( i < n )

    i = i + 1;

    while ( true )
%
%  Select X at random in [-1,+1]
%
      x = rand ( 1, 1 );
%
%  Select Y at random in [0,PDFMAX].
%
      y = pdfmax * rand ( 1, 1 );
%
%  Evaluate Z 
%
      z = 1.0 / sqrt ( pi ) / sqrt ( 1.0 - x .^ 2 ) .^ ( 1.0 / 6.0 );
%
%  Acceptance test:
%
      if ( y <= z )
        break
      end

    end

    s(i) = x;

  end

  return
end

