function prob=bernoull(k,n,p);
% format prob=bernoull(k,n,p);
% Bernoulli trials (binomial pdf):
% prob=probability of k successes in n trials
% with probability p.
% This is also called the binomial distribution
% input:
%     k=scalar or vector, number of successes
%     n=scalar or vector, number of trials
%     p=scalar or vector, probability
%  output:  probability or probability vector
%           rotated to match input size of k.
% reference
% written by E. Gallagher, Environmental Sciences Program
% UMASS/Boston, written: 6/17/93, revised: 6/20/03
% Email: Eugene.Gallagher@umb.edu
[rk,ck]=size(k);[rn,cn]=size(n);[rp,cp]=size(p);
flipk=0;
if rk<ck
   k=k';[rk,ck]=size(k);
   flipk=1;
end
if rn<cn
   n=n';[rn,cn]=size(n);
end
if rn==1 & rk>1
    n=repmat(n,rk,ck);
elseif sum([rn-rk;cn-ck])~=0
      error('Number of elements in k and n must be the same')
end
if rp<cp
   p=p';[rp,cp]=size(p);
end
if rp==1 & rk>1
    p=repmat(p,rk,ck);
elseif sum([rp-rk;cp-ck])~=0
      error('Number of elements in k and n must be the same')
end

o=ones(rk,ck);
lnNf=gammaln(n+o);
lnNmf=gammaln(n-k+o);
lnmf=gammaln(k+o);
if length(p)==1
     prob=exp(lnNf-lnNmf-lnmf+k*log(p)+(n-k)*log(1-p));
else
     prob=exp(lnNf-lnNmf-lnmf+k.*log(p)+(n-k).*log(o-p));
end
if flipk  % return the output vector to match the orientation of the input vector k.
    prob=prob';
end