##############################################################
# R script developped by Christelle Robinet (INRA, Orleans)
# Contact: christelle.robinet@inra.fr
# Date: 30 March 2018
##############################################################



# Function to simulate the dispersal of Monochamus galloprovincialis
#-------------------------------------------------------------------

disp = function(pfmat=0.61,pfimmat=0.45, n=100, agemax=120, nday=120, alpha=2000, reposini=12, repos=1)  	
{
# Inputs:
# pfmat: daily probability of flying for mature beetles
# pfimmat: daily probability of flying for immature beetles
# n: number of insects 
# agemax: adult longevity (in days after adult emergence)
# nday: number of days to simulate the dispersal (if could be < agemax)
# alpha: mean daily dispersal distance (in meters)
# resposini: initial time lag for the dispersal (in days)
# repos: number of days for the rest after one flight (in days)


# tables for the outputs (each column = each day, each line = each beetle):
tabx = matrix(rep(NA,n*(nday+1)),ncol=(nday+1))		# x-coordinates
taby = matrix(rep(NA,n*(nday+1)),ncol=(nday+1))		# y-coordinates
tabx[,1]= 0
taby[,1]= 0
age = 0
compt = 0   		# number of days of rest ; it starts with = 0
exprreposini = F 		# has the initial rest been done ?

for (j in 1:nday)
{
	if ( (age+j-1) < 20 ){nstay = 1-pfimmat; alphaj = (0.67+0.016*(age+j-1))*alpha}	# immatures
	if ( (age+j-1) >= 20 ){nstay = 1-pfmat;  alphaj = alpha}			# matures

	# non flying beetles:
	deadbeetle = ((age+j-1)>agemax)	# insectes morts de "vieillesse"
	exprnoflight = as.logical(rbinom(n=n,size=1,prob=nstay))

if(exprreposini){reposj = repos}; if(!exprreposini){reposj = reposini}; 

if (compt<reposj) {exprrepos = rep(T,n);comptnv = compt+1} 	
# if the rest time is not reached, the insect continues to rest
if (compt==reposj) {exprrepos = rep(F,n);comptnv=0}; 		
# if the rest time has been reached, the insect can fly
if ( (compt==reposj) & (!exprreposini) ) {exprreposini = T}	
# if the initial rest time has been reached, the insect can fly

compt = comptnv

	tabx[exprrepos | exprnoflight | deadbeetle ,(j+1)]= tabx[exprrepos | exprnoflight | deadbeetle ,j]	# non-flying insects stay on the same location
	taby[exprrepos | exprnoflight | deadbeetle ,(j+1)]= taby[exprrepos | exprnoflight | deadbeetle ,j]	
	nflight = sum((!exprrepos)&(!exprnoflight)&(!deadbeetle))
	idflight = (1:n)[(!exprrepos)&(!exprnoflight)&(!deadbeetle)]
	direction = runif(nflight ,min=0,max=2*pi)
	distance = rexp(nflight, rate = 1/alphaj)
	newx = tabx[(!exprnoflight)&(!deadbeetle),j]+distance*cos(direction)
	newy = taby[(!exprnoflight)&(!deadbeetle),j]+distance*sin(direction)
	tabx[idflight,(j+1)]= newx
	taby[idflight,(j+1)]= newy       	
} # end of the loop on j
return(list(tabx=tabx,taby=taby))
} # end of the function
#-------------------------------------------------------------------



# Call of the function for n = 100 insects:
res = disp(n=100) 


# Plot the figure, for instance for day 120 after emergence:
j=120
plot(0,0,pch=19, xlim =c(-50000,+50000),ylim=c(-50000,+50000), xlab="", ylab="", main=paste("Day = ",j), cex.axis=1.3,asp=1)
for (i in 1:100)
{	lines(res$tabx[i,1:j],res$taby[i,1:j],col="darkgreen")  }
points(res$tabx[,j],res$taby[,j],pch=19)