library(spdep)			#Required to use the function cell2nb()

N<-25				#Number of rows and columns in the grid (not the same as N in the paper)
time.steps<-2500		#Number of time steps

amax<-60			#Maximum time steps F. cupressoides can live
bmax<-30			#Maximum time steps Nothofagus spp. can live
disturbed.ratio<-0.5		#Probability a cell is initially disturbed
F.mature<-7			#Age of maturity for F. cupressoides
N.mature<-2			#Age of maturity for Nothofagus spp.
catastrophe.prob<-0.1		#Probability a catastrophe happens in a cell at each time step
neighbor.indices<-cell2nb(N,N,type="rook",torus=TRUE)	#creates a list such that each element
						#stores the indices of the neighbors of that element

F.pop<-as.numeric(runif(N*N)<0.3)	#Is the cell occupied by F. cupressoides?
N.pop<-as.numeric(runif(N*N)<0.5)	#Is the cell occupied by Nothofagus spp.?

s<-as.numeric(runif(N*N)<disturbed.ratio)	#Is the soil in the cell disturbed?
a<-floor((runif(N*N,0,amax+1))*F.pop)		#Occupation time of F. cupressoides in the cell
b<-floor((runif(N*N,0,bmax+1))*N.pop*s)		#Occupation time of Nothofagus spp. in the cell (soil must be disturbed)

timea.colonized<-a			#Creates a matrix of colonization times for F. cupressoides

meanF.density<-array(0,dim=time.steps+1)	#To store density of F. cupressoides on plot
meanN.density<-array(0,dim=time.steps+1)	#To Store density of Nothofagus spp. on plot

meanF.density[1]<-sum(as.numeric(a>0))/(N*N)	#Initial density of F. cupressoides
meanN.density[1]<-sum(as.numeric(b>0))/(N*N)	#Initial density of Nothofagus spp.

for(i in 1:time.steps){
	catastrophe<-as.numeric(runif(length(s))<catastrophe.prob)	#Is there a catastrophe in the cell?
	for(j in 1:length(catastrophe)){
		if(catastrophe[j] == 1){		#If there is a catastrophe...
			s[j]<-a[j]<-b[j]<-0		#Soil is disturbed, all vegetation dies
			timea.colonized<-0		#reset colonization times for these cells
		}
		else{					#If there is not a catastrophe...
			if(s[j] == 0){				#If soil is disturbed...
				b[j]<-0				#Nothofagus spp. cannot grow in the cell
				if(a[j] == 0){			#If cell is not colonized by F. cupressoides then check if it becomes colonized
					if(as.numeric(runif(1)<(sum(as.numeric((a[neighbor.indices[[j]]])>F.mature))/4)) == 1){	
						a[j]<-1
						timea.colonized[j]<-1
					}
					else{ }
				}
				else if(a[j] < amax){		#if cell is colonized by F. cupressoides then add to occupation time
					a[j]<-a[j]+1
					timea.colonized[j]<-timea.colonized[j]+1
				}
				else{				#if occupation time has reach amax then F. cupressoides dies
					a[j]<-0
					timea.colonized[j]<-0
				}
			}
			else{				#If soil is not disturbed...
				if(b[j] == 0){			#If cell is not colonized by Nothofagus spp. then chck if it becomes colonized
					if(as.numeric(runif(1)<(sum(as.numeric((b[neighbor.indices[[j]]])>F.mature))/4)) == 1){
						b[j]<-1
					}
					else {}
				}
				else if(b[j] < bmax){		#If cell is colonized by Nothofagus spp. then add to occupation time
					b[j]<-b[j]+1
				}
				else{				#If occupation time has reach bmax then keep it at bmax
					a[j]<-bmax
				}
				if(a[j] == 0 && b[j] == 0){	#If cell is not colonized by either species check if F. cupressoides colonizes
					if(as.numeric(runif(1)<(sum(as.numeric((a[neighbor.indices[[j]]])>F.mature))/4)) == 1){
						a[j]<-1
						timea.colonized[j]<-1
					}
					else{}
				}				#If cell is colonized by both species then add to occupation time of F. cupressoides
				else if(a[j] < amax && b[j] < bmax){
					a[j] < a[j] + 1
					timea.colonized[j]<-timea.colonized[j]+1
				}
				else if(a[j] == amax || b[j] == bmax){		#If either species has reach its max occupation time then F. cupressoides dies
					a[j]<-0
					timea.colonized[j]<-0
				}
				else{ }
			}
			if(a[j] == 15){		#If F. cupressoides occupation time is 15 then soil in cell is no longer disturbed
				s[j]<-1
			}
			else { }
		}
	}
	meanF.density[i+1]<-sum(as.numeric(a>0))/(N*N)	#Store density of F. cupressoides
	meanN.density[i+1]<-sum(as.numeric(b>0))/(N*N)	#Sotre density of Nothofagus spp.
}

plot(c(0,time.steps),c(0,1),type="n",main="Population Density", xlab="time",ylab="Density")
points(0:time.steps,meanF.density)
points(0:time.steps,meanN.density,pch=2)
legend(1700,0.3,legend=c("F. cupressoides", "Nothofagus spp."), pch=c(1,2))