function out=GenPoisson(in)

%Selects a number from the Poisson Distrubution with mean value of in.
%This is straight from the Wikipedia article: 
%http://en.wikipedia.org/wiki/Poisson_distribution#Generating_Poisson-distributed_random_variables
%
%Example:
%
%   lambda=5;out=GenPoisson(ones([100])*lambda);
%   [counts k]=hist(out(out==out),(0:lambda*3));
%   figure; bar(k,counts/sum(counts));
%   hold on; plot(k,exp(-lambda)*lambda.^k./factorial(k),'r')

if min(in(in==in))<0
    error('values must be zero or positive')
end

out=zeros(size(in));

for ii=1:prod(size(in))
    k=0;
    p=1;
    L=exp(-in(ii));
    while p>L
        k=k+1;
        p=p*rand;
    end
    out(ii)=k-1;
end
