function b = poissd(n,a,t)
%POISSD calculates the value of the Poisson distribution
%
%Call     b = poissd(n,a,t) where
%      n = the maximum integer at which the output is required
%      a = mean events per unit time
%      t = interval over which events are to be counted
%      b = the distribution values for k = 0..n
%
%Note that the distribution is defined for integer arguments only
%
nn = [1,cumprod(1:n)]'; %k!
k = (0:n)';
at = a*t;
eat = exp(-at);
b = eat*at.^k./nn;