function pdf = chebyshev13_pdf ( x )

%*****************************************************************************80
%
%% chebyshev13_pdf() evaluates the Chebyshev1 cube root PDF.
%
%  Formula:
%
%    PDF(X) = ( 1 - x^2 )^(-1/6) * gamm(4/3) / sqrt(pi) / gamma(5/6)
%    for -1 < X < 1
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 August 2016
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, the argument of the PDF.
%
%  Output:
%
%    real PDF, the value of the PDF.
%
  c = gamma ( 4 / 3 ) / sqrt ( pi ) / gamma ( 5 / 6 );

  pdf = zeros ( size ( x ) );

  i = find ( x <= -1.0 | 1.0 <= x );
  j = find ( -1.0 < x & x < 1.0 );

  pdf(i) = 0.0;
  pdf(j) = c * ( 1.0 - x(j) .^ 2 ) .^ ( - 1.0 / 6.0 );

  return
end
