# Growth Models in Practice: Modules 2.2 and 2.3

# Problem 1: trout in a lake with fishing

trout_end_time <- 50
trout_delta_t <- 1 / 12
trout_time <- seq(0, trout_end_time, by = trout_delta_t)

trout_growth_rate <- 0.15
trout_fishing_rate <- 0.08
trout_capacity <- 5000

trout_population <- numeric(length(trout_time))
trout_population[1] <- 400

for (i in 2:length(trout_time)) {
  population_old <- trout_population[i - 1]

  # The logistic term is the net natural change from growth and crowding.
  constrained_growth <- trout_growth_rate * population_old *
    (1 - population_old / trout_capacity)

  # Fishing is an additional outflow.
  caught <- trout_fishing_rate * population_old

  net_rate <- constrained_growth - caught
  trout_population[i] <- population_old + net_rate * trout_delta_t
}

plot(
  trout_time,
  trout_population,
  type = "l",
  lwd = 3,
  col = "darkgreen",
  ylim = c(0, trout_capacity),
  xlab = "Time (years)",
  ylab = "Trout population",
  main = "Trout population with fishing"
)
abline(h = trout_capacity, lty = 2, col = "red")

trout_results <- data.frame(
  time = trout_time,
  population = trout_population
)

reach_2000 <- which(trout_results$population >= 2000)[1]
trout_results[c(reach_2000 - 1, reach_2000), ]

# Problem 2: unconstrained bacterial growth

bacteria_delta_t <- 0.25
bacteria_time <- seq(0, 12, by = bacteria_delta_t)
bacteria_rate <- 0.18

bacteria <- numeric(length(bacteria_time))
bacteria[1] <- 120

for (i in 2:length(bacteria_time)) {
  bacteria_old <- bacteria[i - 1]
  growth <- bacteria_rate * bacteria_old
  bacteria[i] <- bacteria_old + growth * bacteria_delta_t
}

bacteria_exact <- 120 * exp(bacteria_rate * bacteria_time)
bacteria_error <- bacteria - bacteria_exact

bacteria_results <- data.frame(
  time = bacteria_time,
  simulated = bacteria,
  exact = bacteria_exact,
  error = bacteria_error
)

tail(bacteria_results)

reach_500 <- which(bacteria_results$simulated >= 500)[1]
bacteria_results[c(reach_500 - 1, reach_500), ]

plot(
  bacteria_time,
  bacteria,
  type = "l",
  lwd = 3,
  col = "red",
  xlab = "Time (hours)",
  ylab = "Bacteria",
  main = "Unconstrained bacterial growth"
)
lines(bacteria_time, bacteria_exact, lty = 2, lwd = 2)

# Problem 3: caffeine decay

caffeine_half_life <- 5
caffeine_decay_constant <- log(2) / caffeine_half_life

caffeine_delta_t <- 0.25
caffeine_time <- seq(0, 24, by = caffeine_delta_t)

caffeine <- numeric(length(caffeine_time))
caffeine[1] <- 200

for (i in 2:length(caffeine_time)) {
  caffeine_old <- caffeine[i - 1]
  eliminated <- caffeine_decay_constant * caffeine_old
  caffeine[i] <- caffeine_old - eliminated * caffeine_delta_t
}

caffeine_results <- data.frame(
  time = caffeine_time,
  caffeine = caffeine
)

below_50 <- which(caffeine_results$caffeine <= 50)[1]
caffeine_results[c(below_50 - 1, below_50), ]

plot(
  caffeine_time,
  caffeine,
  type = "l",
  lwd = 3,
  col = "purple",
  xlab = "Time (hours)",
  ylab = "Caffeine (mg)",
  main = "Unconstrained decay"
)
abline(h = 50, lty = 2, col = "red")

# Problem 4: constrained recovery of forest biomass

forest_delta_t <- 0.25
forest_time <- seq(0, 80, by = forest_delta_t)
forest_growth_rate <- 0.12
forest_capacity <- 1200

forest_biomass <- numeric(length(forest_time))
forest_biomass[1] <- 80

for (i in 2:length(forest_time)) {
  biomass_old <- forest_biomass[i - 1]
  constraint <- 1 - biomass_old / forest_capacity
  growth <- forest_growth_rate * biomass_old * constraint

  forest_biomass[i] <- biomass_old + growth * forest_delta_t
}

forest_results <- data.frame(
  time = forest_time,
  biomass = forest_biomass
)

forest_results$growth <- forest_growth_rate *
  forest_results$biomass *
  (1 - forest_results$biomass / forest_capacity)

head(forest_results)

plot(
  forest_time,
  forest_biomass,
  type = "l",
  lwd = 3,
  col = "darkgreen",
  ylim = c(0, forest_capacity * 1.05),
  xlab = "Time (years)",
  ylab = "Biomass index",
  main = "Recovery of forest biomass"
)
abline(h = forest_capacity, lty = 2, col = "red")

# When does the forest first reach 75% of carrying capacity?
forest_target <- 0.75 * forest_capacity
reach_forest_target <- which(
  forest_results$biomass >= forest_target
)[1]

forest_results[c(reach_forest_target - 1, reach_forest_target), ]

# When is growth fastest?
fastest_forest_growth <- which(
  forest_results$growth == max(forest_results$growth)
)[1]

forest_results[fastest_forest_growth, ]