function x = gamma_inv( P, m, c )
% x = gamma_inv( P, m, c )
% compute the gamma-distributed variates, x 
% (with mean m and coefficient of variation c) 
% corresponding to non-exceedance probability P :  P = Prob[X <= x] 

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

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

 myeps = 1e-12;
 MaxIter = 100;
 x_old = sqrt(myeps)*ones(size(P));

 for iter = 1 : MaxIter			% iterative Newton's method
	h = (gamma_cdf(x_old, [m,c] ) - P ) ./ gamma_pdf(x_old,m,c);
	x_new = x_old - h;
 	idx = find ( x_new <= myeps );
 	if any (idx)
 	 	x_new(idx) = x_old(idx)/10;
 		h = x_old - x_new;
 	end	
	if max(abs(h)) < sqrt(myeps)	% check convergence
		break;
	end
	x_old = x_new;
 end

 iter

 x = x_new;

% gamma_inv ---------------------------------------------------- 3 March 2015
