# This program simulates the growth of the boxbug population
#   using a stage-structured model.

# DATA

f <- 2.87    # average number of offspring per adult per day
a <- .5   # fraction of adults that survive each day
p <- .5   # fraction of larvae that pupate each day
s <- .25  # fraction of larvae that survive without pupating each day
L0 <- 0   # initial number of larvae
P0 <- 0   # initial number of pupae
A0 <- 2   # initial number of adults
T <- 20   # total number of days

# INITIALIZATION

X <- matrix(0,nrow=3,ncol=T+1)	# X[i,j] is stage i at time j-1
R <- matrix(0,nrow=3,ncol=T)		# R[i,j] is stage i at time j / stage i at time j-1
X[1,1] <- L0
X[2,1] <- P0
X[3,1] <- A0
M <- matrix(0,nrow=3,ncol=3)
M[1,1] <- s
M[2,1] <- p
M[3,2] <- 1
M[3,3] <- a
M[1,3] <- f

# COMPUTATION

for (t in 1:T){	# new time value
	n <- t	# column of the old population values
	k <- t+1	# column of the new population values
	r <- t	# column of the new ratio value
	X[,k] <- M %*% X[,n]
	if (min(X[,n])>0) R[,r] <- X[,k]/X[,n]
		else tr <- t+1	# tr is first column of R that is usable
	}

result <- eigen(M)
lambda1 <- Re(result$values[1])
v1 <- Re(result$vectors[,1]/result$vectors[3,1])

# OUTPUT

par(mfrow=c(1,2))

# plot 1

plot(0:T, X[1,], pch=15, col="green3", xlab="time", ylab="populations")
points(0:T, X[2,], pch=16, col="purple3")
points(0:T, X[3,], pch=17, col="red3")
legend(.05*T,.9*max(X),legend=c("Larvae","Pupae","Adults"),pch=c(15,16,17),col=c("green3","purple3","red3"))
#	[* The legend location may need to be adjusted. *]

# plot 2

plot(tr:T, R[1,tr:T], pch=15, xlim=(c(0,T)), ylim=(c(0,max(R))), col="green3", xlab="time", ylab="growth rates")
points(tr:T, R[2,tr:T], pch=16, col="purple3")
points(tr:T, R[3,tr:T], pch=17, col="red3")
abline(h=lambda1, lty=3)
if (lambda1>0.5*max(R)) ly <- .3 else ly <- .9
legend(.4+.72*T,ly*max(R),legend=c("Larvae","Pupae","Adults"),pch=c(15,16,17),col=c("green3","purple3","red3"))
#	[* The legend location may need to be adjusted. *]

lambda1
v1