function f = gamma_pdf( x, m, c )
% f = gamma_pdf( x, m, c )
% compute the PDF of the gamma probability distribution function
% in terms of it's mean (m) and coefficient of variation (c)


 if (c < 0.5)
   printf('c = %f, lower limit for coefficient of variation is 0.5 \n', c );
 end

 k     = 1/c^2;            % shape parameter
 theta = c^2*m;            % scale parameter

 f = theta^(-k) * x.^(k-1) .* exp(-x/theta) ./ gamma(k);

 f(find(x<0)) = 1e-12; % 0;

% gamma_pdf ---------------------------------------------------- 3 March 2015

