# Module 4.2: Predator–Prey Models
# Squirrels and hawks, adapted from the earlier course example.

# Time and model parameters
delta_t <- 0.01                # months per step
end_time <- 12                 # months
time <- seq(0, end_time, by = delta_t)

squirrel_birth_rate <- 2       # per month
predation_effect <- 0.02       # effect of each hawk on squirrels
hawk_growth_effect <- 0.01     # effect of each squirrel on hawks
hawk_death_rate <- 1.06        # per month

# Stock variables and initial conditions
squirrels <- numeric(length(time))
hawks <- numeric(length(time))
squirrels[1] <- 100
hawks[1] <- 15

# Simulate one time step at a time
for (i in 2:length(time)) {
  squirrels_old <- squirrels[i - 1]
  hawks_old <- hawks[i - 1]

  squirrel_births <- squirrel_birth_rate * squirrels_old
  squirrels_eaten <- predation_effect * squirrels_old * hawks_old
  hawk_growth <- hawk_growth_effect * squirrels_old * hawks_old
  hawk_deaths <- hawk_death_rate * hawks_old

  squirrels[i] <- squirrels_old +
    (squirrel_births - squirrels_eaten) * delta_t
  hawks[i] <- hawks_old +
    (hawk_growth - hawk_deaths) * delta_t
}

# Verify and inspect
results <- data.frame(time, squirrels, hawks)
head(results)
tail(results)

# Populations over time
plot(
  time,
  squirrels,
  type = "l",
  lwd = 3,
  col = "darkgreen",
  ylim = c(0, max(squirrels, hawks)),
  xlab = "Time (months)",
  ylab = "Population",
  main = "Squirrels and hawks over time"
)
lines(time, hawks, lwd = 3, lty = 2, col = "red")
legend(
  "topright",
  legend = c("Squirrels", "Hawks"),
  col = c("darkgreen", "red"),
  lty = c(1, 2),
  lwd = c(3, 3),
  bty = "n"
)

# Largest peaks in this simulated interval
squirrel_peak_position <- which.max(squirrels)
hawk_peak_position <- which.max(hawks)
results[squirrel_peak_position, ]
results[hawk_peak_position, ]

# Predator population versus prey population
plot(
  squirrels,
  hawks,
  type = "l",
  lwd = 2,
  col = "purple",
  xlab = "Squirrels",
  ylab = "Hawks",
  main = "A predator-prey trajectory"
)
points(squirrels[1], hawks[1], pch = 19, col = "red")

# Positive coexistence equilibrium
equilibrium_squirrels <- hawk_death_rate / hawk_growth_effect
equilibrium_hawks <- squirrel_birth_rate / predation_effect
equilibrium_squirrels
equilibrium_hawks