function s = plums ( n )

%*****************************************************************************80
%
%% plums(): average distance to surface of closest of n plums in a spherical pudding.
%
%  Discussion:
%
%    The correct expected value is r/(3*n+1) where r is the radius of the
%    pudding.
%
%  Modified:
%
%    08 May 2019
%
%  Author:
%
%    Paul Nahin
%
%  Reference:
%
%    Paul Nahin,
%    Will You Be Alive 10 Years From Now?,
%    Princeton, 2014,
%    ISBN: 978-0691156804,
%    LC: QA273.25.N344
%
%  Output:
%
%    real s, an estimate for the expected length.
%
  trial_num = 1000000;

  total = 0;
  for trials = 1 : trial_num

    points = zeros ( 1, n );

    for loop2 = 1 : n
%
%  Keep picking random points in the unit cube
%  until you get one that's also in the unit sphere.
%
      while ( true )
        x = 2.0 * rand ( ) - 1.0;
        y = 2.0 * rand ( ) - 1.0;
        z = 2.0 * rand ( ) - 1.0;
        d2 = x^2 + y^2 + z^2;
        if ( d2 < 1.0 )
          break;
        end
      end
%
%  Remember the distance of this point to the center of the sphere.
%
      points(loop2) = sqrt ( d2 );

    end
%
%  Find plum furthest from center (and closest to surface)
%
    total = total + max ( points );

  end

  s = 1.0 - ( total / trial_num );

  return
end
