function rectangle_distance_pdf ( a, b )

%*****************************************************************************80
%
%% rectangle_distance_pdf() plots the PDF for the rectangle distance problem.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 May 2023
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Brahim Gaboune, Gilbert Laporte, Francois Soumis,
%    Expected Distances between Two Uniformly Distributed Random Points
%    in Rectangles and Rectangular Parallelipipeds,
%    The Journal of the Operational Research Society,
%    Volume 44, Number 5, May 1993, pages 513-519.
%    [This reference gives an INCOMPLETE version of PDF!]
%
%    A Mathai, P Moschopoulos, G Pederzoli,
%    Random points associated with rectangles,
%    Rendiconti del Circolo Matematico di Palermo,
%    Volume 47, Number II, pages 163-190, 1999.
%
%  Input:
%
%    real a, b: the length of the two sides.
%
  d = linspace ( 0.0, sqrt ( a^2 + b^2 ), 101 );
%
%  Guarantee a <= b.
%
  if ( b < a )
    t = a;
    a = b;
    b = a;
  end

  s = d .* d;

  term1 = ...
    - 2.0 * d / a^2 / b - 2.0 * d / a / b^2 + pi / a / b + d.^2 / a^2 / b^2;

  term2 = ...
    - 2.0 * d / a^2 / b ...
    - 1.0 / b^2 + 2.0 * asin ( a ./ d ) / a / b ...
    + 2.0 * sqrt ( d.^2 - a^2 ) / a^2 / b;

  term3 = ...
    - 1.0 / b^2 + 2.0 * asin ( a ./ d ) / a / b ...
    + 2.0 * sqrt ( d.^2 - a^2 ) / a^2 / b ...
    - 1.0 / a^2 + 2.0 * asin ( b ./ d ) / a / b ...
    + 2.0 * sqrt ( d.^2 - b^2 ) / a / b^2 ...
    - pi / a / b - d.^2 / a^2 / b^2;

  phi =  term1 .* ( d <= a ) ...
       + term2 .*      ( a < d ) .* ( d <= b ) ...
       + term3 .*                        ( b < d );

  pdf = 2.0 * d .* phi;

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

  return
end

