function x = mypoissrnd(lambda,r,c)
%    x = mypoissrnd(lambda,r,c)
% 
%    Return an R by C matrix of random samples from the Poisson
%    distribution with parameter LAMBDA, which must be a scalar or of
%    size R by C.
%
%    This function works best if LAMBDA is "small" ... ( less than 15 )
%
%    Example: 
%    To view the cumulative sum of a single poisson random variable of 
%    rate paramter 1 use the command ...
%      n=1000; plot(1:n,cumsum(mypoissrnd(1,1,n)), 1:n, 1:n)
%
%    Reference:   http://en.wikipedia.org/wiki/Poisson_distribution

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

 L = exp(-lambda);

 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;

% --------------------------------------------------------------  MYPOISSRND