# preliminaries

Z<-30
H<- matrix (NA, nrow=Z, ncol=Z) 
P<- matrix (NA, nrow=Z, ncol=Z) 
G<- 100 
Beddington<-function(fNt,fPt,fr,falpha,fK,fc){
NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))
PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))
answer<-c(NtPlus1,PtPlus1)
answer
}

random_asign_to_neighbours<-function(to_Add_f, Number_moving, neighbours){
cell<-rep(NA,8)
breaks<- sort(sample(1: Number_moving, 7, replace=TRUE))
for(i in 2:7) {cell[i] <-breaks[i]- breaks[i-1]}
cell[1]<-breaks[1]
cell[8]<- Number_moving - breaks[7]
for(x in 1:8){
to_Add_f[neighbours[x,1],neighbours [x,2]]<- to_Add_f[neighbours[x,1],neighbours [x,2]] + cell[x]
} 
to_Add_f
} # random_asign_to_neighbours function

H[1:Z,1:Z]<-100
P[1:Z,1:Z]<-10
alpha<-0.002
r<-2
c<-1
K<-1000 
proportion_to_migrate<-0.05
neighbours<-matrix(NA, nrow=8, ncol=2)

count<-0

VariancesH<-NULL    #### new
VariancesP<-NULL    #### new


for(g in 1:G) { 
	for(ba in 1:Z){ 
		for(bb in 1:Z){
			answer<-Beddington(H[ba,bb],P[ba,bb],r,alpha,K,c)
			H[ba,bb]<-round(answer[1])
			P[ba,bb]<-round(answer[2])
		} # for bb loop
	} # for ba loop


tomove_Host<-round(H* proportion_to_migrate)
tomove_Parasitoid<-round(P* proportion_to_migrate)
H<- round(H-tomove_Host) 
P<- round(P-tomove_Parasitoid) 
to_add_Host<-matrix(0, nrow=Z, ncol=Z)
to_add_Parasitoid<-matrix(0, nrow=Z, ncol=Z)

for(ba in 1:Z){
	for(bb in 1:Z){  
		neighbours[1,1:2]<- c(ba-1, bb-1) 
		neighbours[2,1:2]<- c(ba-1, bb)
		neighbours[3,1:2]<- c(ba-1, bb+1)
		neighbours[4,1:2]<- c(ba, bb-1)
		neighbours[5,1:2]<- c(ba, bb+1)
		neighbours[6,1:2]<- c(ba+1,bb-1)
		neighbours[7,1:2]<- c(ba+1,bb)
		neighbours[8,1:2]<- c(ba+1,bb+1)
		neighbours[which(neighbours==0)]<- Z 
		neighbours[which(neighbours==Z+1)]<-1 
to_add_Host <- random_asign_to_neighbours (to_add_Host, H[ba,bb], neighbours)
to_add_Parasitoid <- random_asign_to_neighbours (to_add_Parasitoid, P[ba,bb], neighbours) 
	} #end bb loop
} # end ba loop


H<-H + to_add_Host
P<-P + to_add_Parasitoid
VariancesH<-c(VariancesH,var(unlist(as.list(H))))
VariancesP<-c(VariancesP,var(unlist(as.list(P))))
}   # for g loop

#######  new
par(mfrow=c(1,2))
plot(1:G,VariancesH,xlab="Generation",ylab="Variance in hosts per cell",pch=19,col="darkgreen")
plot(1:G,VariancesP,xlab="Generation",ylab="Variance in parasitoids per cell",pch=19,col="darkred")