Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

planck_pdf(x,A,B);
```% PLANCK_PDF evaluates the Planck PDF.  Vector Form!
%
%  Discussion:
%
%    The Planck PDF describes the distribution of frequencies in
%    blackbody radiation at a given temperature T.
%
%    A generalization of the Planck PDF is:
%
%      PDF(X)(A,B) = B**( A + 1 ) * X**A /
%        ( Gamma ( A + 1 ) * Zeta ( A + 1 ) * ( EXP ( B * X ) - 1 ) )
%
%    The standard PDF has A = 3, B = 1.
%
%  Formula:
%
%    PDF(X) = ( 15 / PI**4 ) * X**3 / ( EXP ( X ) - 1 )
%
%  Modified:
%
%    05 March 1999
%
%  Author:
%
%    John Burkardt
%
%  Translated/Modified:
%
%    Jim Huntley,  03/11/04
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    0.0 <= X
%
%    Output, real PDF, the value of the PDF.
%

function [pdf] = planck_pdf(x,A,B);

%persistent gamlnap1 lBap1 lzap1

%if(isempty(lzap1))
lBap1 = log(B^(A+1));
lzap1 = log(zeta(A+1));
gamlnap1 = gammaln(A+1);
%end

sx = size(x,2);

for jx = 1:sx
if (x(jx) < 0)
pdf(jx) = 0;
else
%pdf(jx) = (15.0 / pi^4) * (x(jx))^3 / (exp(x(jx)+eps) - 1);
pdf(jx) = exp(lBap1 + log((x(jx)+eps)^A) - (gamlnap1 + lzap1 + log(exp(B*x(jx)+eps) - 1)));
end
end

return
```