# Module 2.3: Constrained Growth
#
# This example uses Euler's method to simulate a population whose growth is
# limited by a carrying capacity. Time is measured in years.

# 1. Time, parameters, and initial condition

start_time <- 0             # years
end_time <- 100             # years
delta_t <- 1                # years per simulation step

growth_rate <- 0.08         # per year
carrying_capacity <- 500    # individuals
initial_population <- 20    # individuals

# 2. Simulated time points

time <- seq(from = start_time, to = end_time, by = delta_t)

# 3. Stock variable

population <- numeric(length(time))
population[1] <- initial_population

# 4. Euler simulation

for (i in 2:length(time)) {
  population_old <- population[i - 1]
  constraint <- 1 - population_old / carrying_capacity
  growth <- growth_rate * population_old * constraint
  change <- growth * delta_t

  population[i] <- population_old + change
}

# 5. Organize and verify the results

results <- data.frame(time, population)
results$constraint <- 1 - results$population / carrying_capacity
results$growth <- growth_rate * results$population * results$constraint

head(results)

# The first updated population should be 21.536.
results[1:2, ]

# 6. Graph the population and carrying capacity

plot(
  results$time,
  results$population,
  type = "l",
  lwd = 3,
  col = "blue",
  ylim = c(0, carrying_capacity * 1.1),
  xlab = "Time (years)",
  ylab = "Population (individuals)",
  main = "Constrained population growth"
)

abline(
  h = carrying_capacity,
  lty = 2,
  lwd = 2,
  col = "red"
)

legend(
  "bottomright",
  legend = c("Population", "Carrying capacity"),
  col = c("blue", "red"),
  lty = c(1, 2),
  lwd = c(3, 2),
  bty = "n"
)

# 7. When does the population first reach 90% of capacity?

target_population <- 0.90 * carrying_capacity
target_positions <- which(results$population >= target_population)
first_target <- target_positions[1]

results[c(first_target - 1, first_target), ]

# 8. When is growth fastest?

fastest_position <- which(results$growth == max(results$growth))[1]
results[fastest_position, ]