function z = fermi_dirac_sample ( u, v )

%*****************************************************************************80
%
%% fermi_dirac_sample() samples a (continuous) Fermi-Dirac distribution.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 April 2016
%
%  Author:
%
%    Original BASIC version by Frederick Ruckdeschel.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Frederick Ruckdeschel,
%    BASIC Scientific Subroutines,
%    Volume I,
%    McGraw Hill, 1980,
%    ISBN: 0-07-054201-5,
%    LC: QA76.95.R82.
%
%  Input:
%
%    real U, V, the parameters of the distribution.
%    The value of U represents the halfway point for the distribution.
%    Half the probability is to the left, and half to the right, of
%    the value U.  The value of V controls the shape of the distribution.
%    The ratio U/V determines the relative shape of the distribution.
%    Values of U/V in excess of 100 will risk overflow.
%
%  Output:
%
%    real Z, a sample from the Fermi-Dirac distribution.
%    Output values will be nonnegative, and roughly half of them should
%    be less than or equal to U.
%
  iter_max = 1000;

  x = rand ( 1, 1 );
  y = 1.0;
  a = exp ( 4.0 * u / v );
  b = ( x - 1.0 ) * log ( 1.0 + a );

  iter_num = 0;

  while ( true )

    y1 = b + log ( a + exp ( y ) );

    if ( abs ( y - y1 ) < 0.001 )
      break
    end

    y = y1;

    iter_num = iter_num + 1;

    if ( iter_max < iter_num )
      break;
    end

  end

  z = v * y1 / 4.0;

  return
end
