function [HP,JP]=shannon(M,base)
% Calculate Shannon-Wiener diversity H' in (base 2, base 10 ln)
% usage:  function [HP,JP]=shannon(M,base)
%   M is a Station x species data matrix (counts or proportions)
%   base==1 means ln will be used
%   base==2 log-base 2 used
%   base==10 log-base 10 used.
%   HP=Shannon-Wiener's H'
%   JP=Evenness or J'
%   returns NaN if diversity or evenness undefined.
% written by E. Gallagher, Environmental Sciences Program 
% UMASS/Boston, Internet: Gallagher@umbsky.cc.umb.edu
% last revised: 4/22/96
% This will divide each element by the row sum
if nargin==1
  base=1;  % natural logs used if only 1 argument
end
[r,c]=size(M);
M=diag(sum(M'))\M;
z=zeros(r,c);
SPP=sum(~(~M)')';  %SPP(STA,1)
temp=NaN;
HP=temp(ones(r,1),:);
JP=HP;
i=find(M>0);
j=find(SPP>1);
if base == 10
   z(i)=log10(M(i)).*M(i);
   HP=-sum(z');
   JP(j)=(HP(j)./(log10(SPP(j))'))';
elseif base == 2
   z(i)=log2(M(i)).*M(i);
   HP=-sum(z');
   JP(j)=(HP(j)./(log2(SPP(j))'))';
else
   z(i)=log(M(i)).*M(i);
   HP=-sum(z');
   JP(j)=(HP(j)./(log(SPP(j))'))';
end
JP=JP';
