function x = mylognrnd(mu,sigma,r,c)
%    x = mylognrnd(mu,sigma,r,c)
%
%    Return an R by C matrix of random samples from the lognormal
%    distribution with mean MU and standard deviation SIGMA.
%    Both MU and SIGMA must be scalar or of size R by C.
% 
%    If X is a lognormal random variable, then the log of X is normally 
%    distributed with a mean of  log(MU) - 0.5 log (Vx^2+1), 
%    and a variance of  log(Vx^2+1)  where  Vx = SIGMA/MU
%
%    Example: 
%    To view the PDF of a single lognormal random variable of 
%    mean 100 and standard deviation 20, use the command ...
%      n=10000; x=mylognrnd(100,20,1,n); hist(x,n/100); avg_x = sum(x)/n
%
%    Reference: http://en.wikipedia.org/wiki/Log-normal_distribution

 if nargin < 3, [r,c] = size(mu); end
 if any(size(mu) - [r,c] ) & length(mu) ~= 1
    size_of_mu = size(mu)
    r
    c
    error(' mylognrnd: MU must be a scalar or of size R by C')
 end
 if any(size(sigma) - [r,c] ) & length(sigma) ~= 1
    size_of_sigma = size(sigma)
    r
    c
    error(' mylognrnd: SIGMA must be a scalar or of size R by C')
 end
 if length(mu) == 1
    mu    = mu    * ones(r,c);
    sigma = sigma * ones(r,c);
 end

 Vx = sigma./mu;      % coefficient of variation of X

 slx2 = log(Vx.^2+1); % standard deviation of log(X), squared

 z = randn(r,c);      % normally distributed random variable, 0 mean, 1 std-dev

 x = mu .* exp ( -0.5*slx2 + z.*(sqrt(slx2)));

% ----------------------------------------------------------------- MYLOGNRND