neutral = function(maxspec,nocom,maxind,dist,growthfr,pointrate,type_dist_ini,par_dist_ini,pool_ini, nostep, iffig, animation) {
# This function simulates neutral 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: nostep x maxspec matrix with abundance distributions throughout the simulation
# outfreq: nostep x maxspec matrix with species frequencies in communities throughout the simulation
# outsimp: vector with time course of Simpson index of the metacommunity

#   maxspec       - technical limit (slows down calculation)
#   nocom         - number of communities
#   maxind        - size of the community
#   dist          - mortality rate per step
#   growthfr      - proportion of individuals replaced by growth of species present in the community
#                    (migration is 1 - growthfr)      
#   pointrate     - point mutation rate (per individual and step)

#initialization 
#   type_dist_ini - abundance distribution for initialization 
#                     1 geometrical, 2 lognormal, 3 Tokeshi power fraction, 4 random (broken stick) 
#   par_dist_ini  - parameter of the abundance distribution (for geometric and lognormal)
#   pool_ini      - number of species in the initial species pool 

#   nostep        - number of steps

#iffig            - whether to draw figure
#animation        - animation in seconds

  Simpson_multiplier = 5

  noind=matrix(rep(0,nocom*maxspec),nrow=nocom) #communities are rows, species are columns
  names=matrix(rep(0,nocom*maxspec),nrow=nocom)


  for (i in 1:nocom) {
    xx= comm(maxind,type_dist_ini,par_dist_ini,pool_ini,0 );
    noind[i,1:length(xx)] = xx
    names[i,1:length(xx)] = 1:length(xx)
  }

  outdd=NULL;
  outspec=NULL;
  outfreq=NULL;
  outsimp=NULL;

  # initial statistics
   
  totspec=sum(colSums(noind) > 0)
  outspec = c(outspec,totspec)                         #total number of species
  outdd=rbind(outdd,sort(colSums(noind),decreasing=T)) #abundance distribution
  outfreq=rbind(outfreq,colSums(noind>0))              #species frequencies
  outsimp=rbind(outsimp,simpson(colSums(noind)))       #Simpson


  if (iffig) {
    graphics.off()
    a<-system("wmic desktopmonitor get screenwidth",intern=T) 
    scrw = as.numeric(a[2])
      
    dev.new(2,xpos=10,ypos=10,width=5,height=5)
    plot(NULL,xlim=c(0,nostep),ylim=c(0,maxspec), xlab='Time step',ylab="Diversity")
    legend("topleft",lwd=2,col=c("blue","green"),lty=1,legend=c("Number of species", paste("Simpson diversity *",as.character(Simpson_multiplier))))
    dev.new(3,xpos=10++scrw/2,ypos=10,width=5,height=5)
  }

  for (istep in 1:nostep) {
   
    print(istep)
   
    # mortality

    for (i in 1:nocom) {
      for (s in 1:maxspec) {
         noind[i,s] = rbinom(1,noind[i,s], 1-dist);
      }
    }
    names = names * (noind>0);
   
   
    if (nocom == 1)
      sumind = sum(noind)
    else
      sumind = rowSums(noind)
   
    togrow = ceiling((maxind-sumind) * growthfr)
    tomigrate = maxind-sumind-togrow
   
    migrants=NULL
    for (s in 1:maxspec)
      migrants[s]= sum((names==s) * noind)
   
    # local growth

    for (i in 1:nocom) {
      if (sumind[i]) {
         d=rmultinom(1,size=togrow[i],prob=noind[i,])
         noind[i,] = noind[i,] + d
      }
      else {
         tomigrate[i] = maxind
      }
    }
      
    # migration within the metacommunity
   
    if (sum(migrants) > 0) {
      for (i in 1:nocom) {
         d=rmultinom(1,size=tomigrate[i],prob=migrants)
         noind[i,] = noind[i,] + d
         names[i,] = ifelse(noind[i,] == 0 & d > 0, 1:maxspec,names[i,])
      }
    }
   
    # speciation
   
    if (nocom > 1)
      used = as.numeric(colSums(noind) > 0)
    else
      used = as.numeric(noind > 0)
   
    totspec=sum(used)
    if ((pointrate > 0 ) & totspec >= maxspec) {
      print("Maximum community size attained, no speciation possible")
      return()
    }
   
    if (pointrate > 0) {

      for (i in 1:nocom) {
        for (s in 1:maxspec) {
          events = rbinom(1,noind[i,s],pointrate)
          #print(c(i,s,noind[i,s],events))
          if (events) {
            noind[i,s] = noind[i,s] - events
            if (noind[i,s] == 0) 
              names[i,s] = 0
            for (j in 1:events) { 
              if (sum(used) == maxspec) {
                print("Maximum community size attained, no speciation possible1")
                return()
              }
              index=which(used==0)
              noind[i,index[1]] = 1;
              names[i,index[1]] = index[1];
              used[index[1]] = 1;
              #print(c("Speciation event1",events,index[1]))
            }
          }
        }
      }
    }
      
   
    # statistics
   
    totspec=sum(colSums(noind) > 0)
    outspec = c(outspec,totspec)                         #total number of species
    outdd=rbind(outdd,sort(colSums(noind),decreasing=T)) #abundance distribution
    outfreq=rbind(outfreq,colSums(noind>0))              #species frequencies
    outsimp=rbind(outsimp,simpson(colSums(noind)))       #Simpson 
    
    #plotting

    if (iffig) {
      dev.set(2)
      lines(c(istep,istep+1),outspec[istep:(istep+1)], lwd=2,col="blue",lty=1)
      lines(c(istep,istep+1),outsimp[istep:(istep+1)] * Simpson_multiplier, lwd=2,col="green",lty=1)
      dev.set(3)
      barplot(t(noind),col=rainbow(maxspec),main=paste(istep))
      Sys.sleep(animation)
    }
  }

  outddend = sort(colSums(noind),decreasing=T)
  list(outddend=outddend,outspec=outspec,outdd=outdd,outfreq=outfreq,outsimp=outsimp)
}


simpson = function(x) {
  1/sum((x / sum(x))^2)
}