source("comm.R")
source("distr.R")
source("neutral.R")


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

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

nostep = 2000;       # number of steps
iffig = 1
animation=0.

a=neutral(maxspec,nocom,maxind,dist,growthfr,pointrate,type_dist_ini,par_dist_ini,pool_ini, nostep, iffig, animation) 

plot(a$outspec,type="l",lwd=2,col="blue",xlab="Time step",ylab="Number of species")
plot(a$outsimp,type="l",lwd=2,col="blue",xlab="Time step",ylab="Simpson diversity")
plot(a$outspec,a$outsimp,type="l",lwd=2,col="blue",ylab="Simpson diversity",xlab="Number of species")
plot(a$outddend[a$outddend>0],log="y",type="b",lwd=2,col="blue",xlab="species",ylab="Log abundance")