function x = Poisson_rnd(T,r,c)
%    x = Poisson_rnd(T,r,c)
% 
%    Return an R by C matrix of random samples from the Poisson
%    distribution with return period T, which must be a scalar or of
%    size R by C.
%
%    This function works best if T is "large" ... ( more than than 0.1 )
%
%    Example: 
%    To view the cumulative sum of a single Poisson random variable of 
%    rate paramter 1 use the commands ...
%    n=1000; plot(1:n,cumsum(Poisson_rnd(1,1,n)), 1:n, 1:n)
%
%    Reference:   http://en.wikipedia.org/wiki/Poisson_distribution

% H.P. Gavin, Dept. Civil & Environ. Eng'g, Duke Univ.

 if ( nargin < 2 ), [r,c] = size (T); end
 if ( any(size(T) - [r,c]) ) && ( length(T) ~= 1 )
    size_of_T= size(T)
    r
    c
    error(' Poisson_rnd: T must be a scalar or of size R by C')
 end

 L = exp(-1./T);

 p  = ones(r,c);
 x  = zeros(r,c);
 ij = find(p>=L);

 while length(ij) > 0
   p(ij)  =  p(ij) .* rand(size(p(ij)));
   x(ij)  =  x(ij) + 1;
     ij   =  find(p>=L);
 end
 x = x - 1;

% ----------------------------------------------------------------  Poisson_rnd
% 2019-04-20
