function mu = box_distance_mean ( a, b, c )

%*****************************************************************************80
%
%% box_distance_mean() returns the exact mean of the box distance.
%
%  Discussion:
%
%    The box has sides of length a <= b <= c.
%
%    For a = b = 1, we should get 
%    mu = ( 4 + 17 sqrt(2) - 6 sqrt(3) + 21 log(1+sqrt(2)) + 84 log(1+sqrt(3))
%      - 42 log(2) - 7 pi ) / 105 = 0.661707182...
%
%    There would seem to be an error in my implementation, since
%    box_distance_mean(1,1,1) is coming out to be 0.678662...
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 May 2023
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Johan Philip,
%    The probability distribution of the distance between two random
%    points in a box,
%    https://people.kth.se/~johanph/habc.pdf
%
%  Input:
%
%    real a, b, c: the length of the sides.
%
%  Output:
%
%    real mu: the exact mean distance between two random points in the box.
%
  r =  sqrt ( a^2 + b^2 + c^2 );
  r1 = sqrt (       b^2 + c^2 );
  r2 = sqrt ( a^2       + c^2 );
  r3 = sqrt ( a^2 + b^2       );

  mu = r / 15.0 ...
     - 7.0 * ( r - r1 ) * r1^2 / 90.0 / a^2 ...
     - 7.0 * ( r - r2 ) * r2^2 / 90.0 / b^2 ...
     - 7.0 * ( r - r3 ) * r3^2 / 90.0 / c^2 ...
     + 4.0 * ( a^7 + b^7 + c^7 - r1^7 - r2^7 - r3^7 + r^7 ) ...
     / 315.0 / a^2 / b^2 / c^2 ...
     + ( b^6 * asinh ( a / b ) + c^6 * asinh ( a / c ) ...
     - r1^2 * ( r1^4 - 8.0 * b^2 * c^2 ) * asinh ( a / r1 ) ) ...
     / 30.0 / a / b^2 / c^2 ...
     + ( c^6 * asinh ( b / c ) + a^6 * asinh ( b / a ) ...
     - r2^2 * ( r2^4 - 8.0 * c^2 * a^2 ) * asinh ( b / r2 ) ) ...
     / 30.0 / a^2 / b / c^2 ...
     + ( a^6 * asinh ( c / a ) + b^6 * asinh ( c / b ) ...
     - r3^2 * ( r3^4 - 8.0 * a^2 * b^2 ) * asinh ( c / r3 ) ) ...
     / 30.0 / a^2 / b^2 / c ...
     - 2.0 * ( a^4 * asin ( b * c / r2 / r3 ) ...
             + b^4 * asin ( c * a / r3 / r1 ) ...
             + c^4 * asin ( a * b / r1 / r2 ) ) / 15.0 / a / b / c;

  return
end

