function x = Lapalce_rnd(mX,sX,r,c)
% x = Laplace_rnd(mX,sX,r,c)
%
% generate a matrix of a random variable distribued according to the 
% Laplace distribution with parameters mX and sX
%
%  INPUT         DESCRIPTION
%   mX           mean    of the Laplace distribution          1 x 1   or  r x c
%   sX           std.dev of the Laplace distribution          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 Laplace distribution                       r x c
%
%  EXAMPLE 
%    To plot the PDF of a single Laplace random variable, ... 
%    n=10000; x=Laplace_rnd(1,1,1,n); hist(x,n/100,1); avg_x = sum(x)/n
 
%    References:  http://en.wikipedia.org/wiki/Exponential_distribution
%                 http://www.math.uah.edu/stat/special/Laplace.html

 if any(sX<=0) || any(sX==Inf)
    error(' Laplace_rnd: sX must be greater than zero')
 end
 
 sr2 = sqrt(2);
 
 if nargin == 3
   u = norm_cdf(r);
   [r,c] = size(u);
 end
 if nargin == 4
   u = rand(r,c);
 end 

 if (prod(size(mX)) == 1), mX = mX * ones(r,c); end
 if (prod(size(sX)) == 1), sX = sX * ones(r,c); end

 in = find(u<=0.5);
 ip = find(u>0.5);

 x(in) = mX(in) + sX(in)/sr2 .* log(2*u(in));
 x(ip) = mX(ip) - sX(ip)/sr2 .* log(2*(1-u(ip)));

% x = mX - (sX/sr2) .* u .* log(2.0 - 2.0*abs(u));  % ???

% --------------------------------------------------------------  Laplace_rnd
% HP Gavin 2021-04-26
