function [cnt, pd1, pd2] = poisscnt2(rate1,rate2,interval,maxcount)
% poisscnt2: function which returns the probability density for the 
% distribution of the spike counts of two Poisson process in a fixed time
% interval. If called without output arguments, the probability
% distributions are plotted.
%
%	[cnt1, pd1, pd2] = poisscnt2(rate1,rate2,interval,maxcount)
%
%	where 
%	  rate1, rate2 = mean rate of the Poisson processes [Hz]
%	  interval = the length of the interval (in msec)
%	  maxcount = maximal count taken into account
%
% The return parameters are:
%	cnt = vectors of spike counts.
%	pd1, pd2 = corresponding probability distributions.
%
%

if ( nargin ~= 4 )
  disp(' ');
  disp('usage: poisscnt2(rate1,rate2,interval,maxcount) ');
  disp('       for more information type "help poisscnt" in the main');
  disp('       matlab window');
  disp(' ');
  return;
end;

%initializes the various variables
lambda1 = rate1*interval*1e-3;
lambda2 = rate2*interval*1e-3;

cnt = (0:1:maxcount)';
pd1 = zeros(maxcount+1,1);

for k = 1:maxcount+1
  if ( (prod(1:k-1) ~= Inf) & (lambda1^(k-1) ~= Inf) )
    pd1(k,1) = lambda1^(k-1)*exp(-lambda1)/prod(1:k-1);
  else 
    pd(k,1) = 0;
  end;
end;

pd2 = zeros(maxcount+1,1);

for k = 1:maxcount+1
  if ( (prod(1:k-1) ~= Inf) & (lambda2^(k-1) ~= Inf) )
    pd2(k,1) = lambda2^(k-1)*exp(-lambda2)/prod(1:k-1);
  else 
    pd2(k,1) = 0;
  end;
end;