Fig_04

#---- Nicholson & Bailey
# fig_04

# variables
N <- numeric()
P <- numeric()

# coeficientes
lambda <- 1.2
c <- 4
a <- 0.005
R <- 5

# condiciones simulación
ngen <- 200    # generaciones
N[1] <- 70     # N inicial
P[1] <- 60     # P inicial

# simulación
for (t in 1:ngen) {
  f <- exp(-a * P[t])
  N[t + 1] <- lambda * N[t] * f
  P[t + 1] <- lambda * c * N[t] * (1 - f) + R
}

# equilibrio
Peq <- log(lambda) / a
Neq <- (Peq - R) / (c * (lambda - 1))

# gráfico
pdf("fig_04.pdf", width = 7.5, height = 6)
par(pty = "s",
    bty = "n",
    mar = c(5, 2, 1, 1),
    xaxs = "i",
    yaxs = "i",
    las = 1)
plot(N, P,
    type = "l",
    lty = 2,
    col = "lightblue",
    xlim = c(0, 80),
    ylim = c(0, 120))
points(N, P,
    cex = 0.8,
    pch = 16,
    col = "royalblue")
abline(h = Peq,
    col = "orange",
    lwd = 0.5)
abline(v = Neq,
    col = "orange",
    lwd = 0.5)
abline(h = R,
    col = "royalblue",
    lwd = 0.5)
polygon(c(0, 0, 80, 80),
    c(0, R, R, 0),
    density = 10,
    col = "grey")
axis(1)
axis(2)
dev.off()