function x = triangular_inv(P,a,b,c)
%    x = triangular_inv(P,a,b,c)
% 
%    Return the inverse CDF of the triangular distribution with limits 
%    a < X < b, and mode c, which must be scalars. 
%
%    Reference:   http://en.wikipedia.org/wiki/Triangular_distribution

 P(find(P<=0)) = eps;
 P(find(P>=1)) = 1-eps;

 Fc = (c-a) ./ (b-a);		% CDF at x=c

 x = b - sqrt((1-P).*(b-a).*(b-c));			% ... if 0 < P < Fc

 idx = find(P<Fc);

 x(idx) = a + sqrt(P(idx).*(b-a).*(c-a));		% ... if Fc < P < 1 

% ----------------------------------------------------------  triangular_inv
