function [ z, seed ] = fermi_dirac_sample ( u, v, seed ) %*****************************************************************************80 % %% FERMI_DIRAC_SAMPLE samples a (continuous) Fermi-Dirac distribution. % % Licensing: % % This code is distributed under the GNU LGPL 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. % % Parameters: % % 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. % % Input/output, integer SEED, a seed for the random % number generator. % % 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, seed ] = r8_uniform_01 ( seed ); y = 1.0; a = exp ( 4.0 * u / v ); b = ( x - 1.0 ) * log ( 1.0 + a ); iter_num = 0; while ( 1 ) 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