function cube_distance_pdf ( )

%*****************************************************************************80
%
%% cube_distance_pdf plots the PDF for the cube distance problem.
%
%  Discussion:
%
%    Determining the correct MATLAB expression for this formulation is an
%    excellent test of patience and persistence!
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 December 2018
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eric Weisstein,
%    Cube Line Picking,
%    From MathWorld--A Wolfram Web Resource.,
%    http://mathworld.wolfram.com/CubeLinePicking.html
%
  d = linspace ( 0.0, sqrt ( 3.0 ), 101 );

  pdf = ( 0.0 <= d ) ...
      .* ( d <= 1.0 ) ...
      .* ( - d .^ 2 ) ...
      .* ( ( d - 8.0 ) .* d .^ 2 + pi * ( 6.0 * d - 4.0 ) ) ...
      + ( 1.0 < d ) ...
      .* ( d <= sqrt ( 2.0 ) ) ...
      .* 2.0 .* d .* ...
      ( ...
        ( d .^ 2 - 8.0 * sqrt ( d .^ 2 - 1.0 ) + 3.0 ) .* d .^ 2 ...
        - 4.0 * sqrt ( d .^ 2 - 1.0 ) ...
        + 12.0 * d .^ 2 .* asec ( d ) ...
        + pi * ( 3.0 - 4.0 * d ) - 0.5 ...
      ) ...
      + ( sqrt ( 2.0 ) < d ) ...
      .* ( d <= sqrt ( 3.0 ) ) ...
      .* d .* ...
      ( ...
        ( 1.0 + d .^ 2 ) ...
        .* ( 6.0 * pi + 8.0 * sqrt ( d .^ 2 - 2.0 ) - 5.0 - d .^ 2 ) ...
        - 16.0 .* d .* acsc ( sqrt ( 2 - 2 ./ d .^ 2 ) ) ...
        + 16.0 .* d .* atan ( d .* sqrt ( d .^ 2 - 2.0 ) ) ...
        - 24.0 .* ( d .^ 2 + 1.0 ) .* atan ( sqrt ( d .^ 2 - 2.0 ) ) ...
      );

  plot ( d, pdf, 'r-', 'linewidth', 2 );
  grid ( 'on' );
  xlabel ( '<-- Distance -->' );
  ylabel ( '<-- Probability -->' );
  title ( 'PDF for distance between pairs of random points in unit cube' );

  return
end
