function [HB,E]=brillouin(DATA)
% Calculation of Brillouin species diversity.
% HB=(ln N!-sum(ln ni!))/N
% format [HB,E]=Brillouin(DATA)
% Input: DATA is a sample by species matrix, 
% Output: HB is Brillouin sample diversity
%         E is the Evenness index, Pielou's V
%    Hurlbert SM (1971)  The non-concept of species diversity:  a
%           critique and alternative parameters.  
%          Ecology 52:  577-586.
%    Magurran AE (1988) Ecological diversity and its measurement.
%            Princeton Unversity Press.
% written 4/5/94 by E. Gallagher,
% Dept of Environmental, Earth & Ocean Sciences
% UMASS/Boston, Eugene.Gallagher@umb.edu
% last revised 9/27/2010
[nsamp,nspec]=size(DATA);
N=sum(DATA')';
O=ones(nsamp,nspec);
% find the number of species.
S=sum(~(~DATA)')';
lnNfact=gammaln(N+1);
sumlnnifact=sum(gammaln(DATA+O)')';
HB=(lnNfact-sumlnnifact)./N;
if nargout>1
      NS=floor(N./S);
      r=N-S.*NS;
       HBMAX=(ones(nsamp,1)./N).* ...
               (lnNfact-((S-r).*gammaln(NS+1)+...
      r.*gammaln(NS+2)));
      E=HB./HBMAX;
end