function [outddend,outspec,outdd,outfreq,lifespan,outsimp] = hubbell(iffig,maxspec,nocom,dist,growthfr,outside,maxind,pointrate,splitrate,type_dist_ini,par_dist_ini,pool_ini, nostep )
% This function simulates neutral processes of community dynamics and speciation

%Output (returns these five variables): 
% outddend: vector of abundance distribution at the end of simulation
% outspec: vector with time course of species richness of the metacommunity
% outdd: maxspec x nostep matrix with abundance distributions throughout the simulation
% outfreq: maxspec x nostep matrix with species frequencies in communities throughout the simulation
% lifespan: vector of durations of individual species (in time steps)
% outsimp: vector with time course of Simpson index of the metacommunity

if (exist('maxspec') == 0)

   maxspec = 200; % just a technical limit (slows down calculation)
   nocom=10 ;     % number of communities
   dist=0.3;      % mortality rate per step
   growthfr=0.6;  % proportion of individuals replaced by growth of species present in the community
   outside = 0;   % proportion of immigrants coming from the outside species pool. For Hubbell's model should be == 0
   maxind=100;    % size of the community
   pointrate = 0.005; %point mutation rate (per individual and step)
   splitrate = 0.; %species splitting (allopatric speciation) rate (per species and step)

   %initialization (also used as species pool if outside > 0)

   type_dist_ini=4;  % abundance distribution for initialization and sampling from the species pool
                  % 1 geometrical, 2 lognormal, 3 Tokeshi power fraction, 4 random (broken stick) 
   par_dist_ini=0.9; % parameter of the abundance distribution
   pool_ini = 1;     % number of species in the species pool 

   nostep = 1000;     % number of steps
end


%**************************************************************************

noind=zeros(nocom,maxspec);
names=zeros(nocom,maxspec);


for i=1:nocom
   xx= fliplr(comm(maxind,type_dist_ini,par_dist_ini,1,pool_ini,0 ))';
   noind(i,:) = [xx' zeros(1,maxspec - pool_ini)];
   names(i,:) = (1:maxspec) .* (noind(i,:)>0);
end

outdd=[];
outspec=[];
outfreq=[];
outsimp=[];

size(noind);

if (iffig)
   scrsz = get(0,'ScreenSize');
   figure('Position',[10 scrsz(4)/2-30 scrsz(3)/2-30 scrsz(4)/2-60], 'Name', 'Number of species')
   axis([0, nostep,0, maxspec]);
   hold on;
   fig2=gcf;
   figure('Position',[13*scrsz(4)/20 scrsz(4)/2-30 scrsz(3)/2-30 scrsz(4)/2-60], 'Name', 'Metacommunities')
   bar(noind,'Stacked');
   fig1=gcf;
   figure('Position',[13*scrsz(4)/20 20 scrsz(3)/2-30 scrsz(4)/2-60], 'Name', 'Abundance distribution')
   fig3=gcf;
   hold on;
   figure('Position',[10 20 scrsz(3)/2-30 scrsz(4)/2-60], 'Name', 'Simpson index')
   axis([0, nostep,0, 1]);
   hold on;
   fig4=gcf;
      
end



for istep=1:nostep
   
   istep
   
   % mortality
   for i=1:nocom
      for s=1:maxspec
         noind(i,s) = rbinom(1,noind(i,s), 1-dist);
      end
   end
   names = names .* (noind>0);
   
   %'mort'
   %noind

   
   sumind = sum(transpose(noind));
   togrow = ceil((maxind-sum(transpose(noind))) * growthfr);
   tomigrate = floor((maxind-sum(transpose(noind))) * (1 - growthfr) * (1 - outside));
   tocome = floor((maxind-sum(transpose(noind))) * (1 - growthfr) * (outside));
   
   for s=1:maxspec
      migrants(s)= sum(sum((names==s) .* noind));
   end

   % local growth

   for i=1:nocom
      %[i sumind(i)];
      if (sumind(i))
         d=cumsum(noind(i,:))./sumind(i);
         n = zeros(maxspec,1);
         for s=1:togrow(i)
            r = rand(1);
            for j=1:maxspec
               if r<d(j)
                  n(j) = n(j)+1;
                  break;
               end
            end
         end
         %n2=noind(i,:)
         noind(i,:) = noind(i,:) + n';
      else
         tomigrate(i) = ceil(maxind* (1 - outside));
         tocome(i) = floor(maxind* outside);

      end
   end
   
   %'growth'
   %noind
   
   % migration within the metacommunity
   
   if (outside ~= 1 & sum(migrants)>0)
      d=cumsum(migrants)./sum(migrants);
      for i=1:nocom
         %tomigrate(i); 
         %noi=noind(i,:)
         n = zeros(maxspec,1);
         for s=1:tomigrate(i)
            r = rand(1);
            for j=1:maxspec
               if r<d(j)
                  n(j) = n(j)+1;
                  break;
               end
            end
         end
         %n2=noind(i,:)
         noind(i,:) = noind(i,:) + n';
         %n3=noind(i,:)
         %n=n';
         %s=sum(noind(i,:))
      end
   end
   %'migrat'
   %noind
   
   % immigration from the species pool
   
   if (outside ~= 0)
      for i=1:nocom
         n = [comm(tocome(i),type_dist_ini,par_dist_ini,1,pool_ini,0 )' zeros(1,maxspec - pool_ini)];
         noind(i,:) = noind(i,:) + n;
         %n3=noind(i,:)
         %n=n';
         %s=sum(noind(i,:))
      end
   end
   
   % speciation
   
   totspec=sum(sum(noind(1:nocom,:)>0)>0);
   if ((pointrate > 0 | splitrate > 0) & totspec == maxspec)
      'Maximum community size attained, no speciation possible'
      return;
   end
   
   if (pointrate > 0) %point mutation

      used=sum(noind>0)>0;
      for i=1:nocom
         for s=1:maxspec
            for j=1:noind(i,s) 
               if rand<pointrate
                  index=find(used==0);
                  noind(i,s) = noind(i,s)-1;
                  noind(i,index(1))=1;
                  names(i,index(1))=s;
                  totspec=totspec+1;
                  used(index(1))=1;
                  index(1);
                  'Speciation event1';
                  %totspec
                  %maxspec
                  if (totspec == maxspec)
                     'Maximum community size attained, no speciation possible 1'
                     return;
                  end
               end
            end
         end
      end
   end
   
      
   if (splitrate > 0) %species splitting (allopatric speciation)
      used=sum(noind>0)>0;
      for s=1:maxspec
         %s
         %ss=sum(noind(:,s)>1)
         if rand<splitrate & sum(noind(:,s)>1)>1
            'Speciation event2';
            index=find(used==0);
            index=index(1);
            if length(index) == 0
               'Maximum community size attained, no speciation possible 2'
               return;
            end           
            for i=1:nocom
               if rand < 0.5
                  noind(i,index) = noind(i,s);
                  %noi=noind(i,index) 
                  noind(i,s) = 0;
                  names(i,index) = s;
                  names(i,s) = 0;
               end
            end
         end
      end
   end
   
   % statistics
   
   %noind
   %sum(sum(noind))
   
   if (nocom > 1)
      outdd=[outdd;sum(noind)];                          %abundance distribution
      simptemp=sum( (sum(noind)./sum(sum(noind))).^2);      %Simpson index
      outsimp=[outsimp simptemp];
      totspec=sum(sum(noind(1:nocom,:)>0)>0);
      outspec=[outspec totspec];                         %total number of species
      outfreq=[outfreq;sum(noind>0)./nocom ];            %species frequencies
   else
      outdd=[outdd;noind];                          %abundance distribution
      simptemp=sum( (noind./maxind).^2)             %Simpson index
      outsimp=[outsimp simptemp];
      totspec=sum(noind>0);
      outspec=[outspec totspec];                    %total number of species
      outfreq=[outfreq;sum(noind>0)];               %species frequencies
   end
   
      
   
   
   if (iffig)
      
         
      figure(fig1);
      clf
      bar(noind,'Stacked');
      hold on
      pause(0.0);
   
      figure(fig2);
      if (istep>1)
         plot(istep-1:istep,[outspec(istep-1) outspec(istep)]);
      else
         plot(0:1,[pool_ini outspec(istep)])
      end
      
      figure(fig4);
      if (istep>1)
         plot(istep-1:istep,[outsimp(istep-1) outsimp(istep)]);
      else
         %plot(0:1,[0 outsimp(istep)])
      end
 
      
      
      
      figure(fig3)
      if (istep>1)
         plot(xx);
      end
      if (nocom > 1)
         xx = sort(sum(noind));
      else
         xx = sort(noind);
      end
      
      xx = fliplr(xx);
      
      xx=[log(xx(find(xx))) -0.5* ones(1,maxspec-length(find(xx)))];
      
      plot(xx,'Color','r');
   end

   
end
outddend=sort(sum(noind)); %final abundance distribution
outddend=fliplr(outddend);

%lifespan calculation
lifespan=[];
df=diff(outdd>0);
for s=1:maxspec
   ss=sum(outdd(1:nostep,s));
   if (ss>0)
      df(1:nostep-1,s);
      if (sum(df(1:nostep-1,s)) == 0)
         lifespan = [lifespan nostep];
      else
         x=find(df(1:nostep-1,s)~= 0);
         dif=diff([x' nostep]);
         lifespan = [lifespan dif(1:2:length(dif))];
      end
   end
end
lifespan = fliplr(sort(lifespan));

%'Final abundance distribution:'