function x = corr_logn_rnd( medX , covX , R, N )
% x = corr_logn_rnd( medX , covX , R, N )
%
% generate a sample of m observations of n log normal random variables to the 
% with parameters  medX (n by 1)  and  covX (n by 1)
% the standard normal deviates of X are correlated with matrix R
%
%  INPUT         DESCRIPTION
%   medX         median of the lognormal distribution                     n x 1
%   covX         c.o.v. of the lognormal distribution                     n x 1
%   R            correlation matrix of standard normal deviates           n x n
%   N            number of observations                                   1 x 1
%
%    Values of  medX  and  covX  must be greater than zero.  
%
%  OUTPUT
%   x            sample of m observations of n correlated lognormal r.v's n x m
%
%    If X is a lognormal random variable, then the log of X is normally 
%    distributed with a mean of  log(medX) and a variance of log(1+covX^2)
%
%  EXAMPLE 
%    To plot the PDF of a single lognormal random variable, ... 
%    n=10000; x=logn_rnd(1,0.5,1,n); hist(x,n/100,1); avg_x = sum(x)/n
 
%    Reference: http://en.wikipedia.org/wiki/Log-normal_distribution

 if nargin ~= 4
    help corr_logn_rnd
    return
 end

 medX = medX(:);           % make medX a column vector
 covX = covX(:);           % make covX a column vector

 n = length(medX);         % number of random variables

 if any(diag(R)~=1), error(' corr_logn_rnd: diagonal of R must equal 1'); end
 if any(R>1),  error(' corr_logn_rnd: R values cannot be greater than 1'); end
 if any(R<-1), error(' corr_logn_rnd: R values cannot be less than 1');   end
 if size(medX,1)~=n, error(sprintf(' b must have %d rows, zero columns',n)); end
 if size(covX,1)~=n, error(sprintf(' q must have %d rows, zero columns',n)); end
 if any(medX<0), error(' medX must be positve'); end
 if any(covX<0), error(' covX must be positve'); end

 VlnX = log(1+covX.^2);	   % the variance of log(X)

% eigenvalues and eigenvectors of the correlation matrix R
 [eVec,eVal] = eig(R);
 if any(eVal<0), error(' corr_logn_rnd:  R must be positive definite'); end

% m observations of n uncorrelated standard normal r.v's
 Z = randn(n,N); 

% N observations of n correlated standard normal r.v's
 Y = eVec * sqrt(eVal) * Z;  

% N observations of n correlated lognorm r.v's
 x = exp ( log(medX) + Y.*(sqrt(VlnX)*ones(1,N)));

% ------------------------------------------------------------- corr_logn_rnd
% HP Gavin  2016-04-28, 2021-04-21, 2022-02-12
