function x = triangular_rnd(A,B,C,r,c)
% x = triangular_rnd( A, B, C, r,c)
%
% generate a matrix of a random variable distribued according to the 
% triangular distribution with parameters  A, B, and C
%
%  INPUT         DESCRIPTION
%   A            lower bound of X                             1 x 1   or  r x c
%   B            upper bound of X                             1 x 1   or  r x c
%   C            mode of X                                    1 x 1   or  r x c
%   r            number of rows in x                                      1 x 1
%                sample of standard normal variables, possibly correlated r x c
%   c            number of cols in x                                      1 x 1
%
%  OUTPUT
%   x            sample of the extreme type I distribution                r x c
%
%  EXAMPLE 
%    To plot the PDF of a single extreme type I random variable, ... 
%    n=10000; x=triangular_rnd(1,4,2,1,n); hist(x,n/100,1); avg_x = sum(x)/n

%    Reference:   http://en.wikipedia.org/wiki/Triangular_distribution

 if any(B<A), error(' triangular_rnd:  A must be less than B'); end
 if any(C<A), error(' triangular_rnd:  A must be less than C'); end
 if any(B<C), error(' triangular_rnd:  C must be less than B'); end

 if nargin == 4
   u = norm_cdf(r);
   [r,c] = size(u);
 end
 if nargin == 5
   u = rand(r,c);
 end

 if (prod(size(A)) == 1), A = A * ones(r,c); end
 if (prod(size(B)) == 1), B = B * ones(r,c); end
 if (prod(size(C)) == 1), C = C * ones(r,c); end

 Fc = (C-A) ./ (B-A);		% CDF at x = C

 x = B - sqrt((1-u).*(B-A).*(B-C));			% ... if 0 < u < Fc 

 idx = find(u<Fc);

 x(idx) = A(idx) + sqrt(u(idx).*(B(idx)-A(idx)).*(C(idx)-A(idx)));

% ----------------------------------------------------------  triangular_rnd
% HP Gavin 2016-04-29, 2021-04-26
