####################
# Global Constants #
####################

# age at which birds in captivity begin breeding (become adults)
captive_begin_breeding <- 7
# age at which birds in the wild begin breeding (become adults)
wild_begin_breeding <- 8
# if the captive population exceeds this amount, then the entire year one
# population will be released
max_in_captivity <- 150
# carrying capacity in the wild
wild_carrying_capacity <- 400
# percentage of children which will make it to year 1 in captivity,
# per mating pair
captive_productivity_rate <- 0.9
# percentage of children which will make it to year 1 in captivity,
# per mating pair (with double clutching)
captive_productivity_rate_dc <- 1.8
# percentage of children which will make it to year 1, per mating pair
wild_productivity_rate <- 0.3729

################################
# Global Simulation Parameters #
################################

# number of years to simulate
simulation_years <- 50

######################
# Derived Parameters #
######################

# rate at which captive immature birds will become adults
# an approximation which assumes an equal distribution of immature bird ages
captive_immature_to_adult_rate <- 1 / (captive_begin_breeding - 2)
# rate at which wild immature birds will become adults
# an approximation which assumes an equal distribution of immature bird ages
wild_immature_to_adult_rate <- 1 / (wild_begin_breeding - 2)

# Simulates the Condor population over time using the given parameters.
# * `rd1` - death rate of year one birds
# * `rdI` - death rate of immature birds
# * `rdA` - death rate of adult birds
simulate_condor_population <- function(rd1, rdI, rdA) {
  # captive population groups
  c_Y0 <- c(10)
  c_Y1 <- c(10)
  c_I <- c(10)
  c_A <- c(10)

  # wild population groups
  w_Y0 <- c()
  w_Y1 <- c()
  w_I <- c()
  w_A <- c()

  years <- seq(1, simulation_years)

  for (t in 2:simulation_years) {
    # model captive population changes (year 1, immature, and adult)
    c_Y1[t] <- c_Y0[t - 1] - # last year's population
      rd1 * c_Y0[t - 1] # outflow to deaths

    c_I[t] <- c_I[t - 1] - # last year's population
      rdI * c_I[t - 1] - # outflow to deaths
      captive_immature_to_adult_rate * c_I[t - 1] + # outflow to adult
      c_Y1[t - 1] # inflow from year 1

    c_A[t] <- c_A[t - 1] - # last year's population
      rdA * c_A[t - 1] + # outflow to deaths
      captive_immature_to_adult_rate * c_I[t - 1] # inflow from immature

    # model new births in captive population
    c_Y0[t] <- captive_productivity_rate * (c_A[t] / 2)

    # sum of all population groups
    captive_population <- c_Y0[t] + c_Y1[t] + c_I[t] + c_A[t]
    released_to_wild <- if (captive_population > max_in_captivity) {
      c_Y1[t]
    } else {
      0.8 * c_Y1[t]
    }
    # subtract the released birds from the captive year one population
    c_Y1[t] <- c_Y1[t] - released_to_wild
  }

  captive_population <- c_Y0 + c_Y1 + c_I + c_A
  plot(
    years,
    c_Y0,
    ylim = c(0, 50),
    col = "red",
    main = "Condor Population in Captivity",
    xlab = "Time (years)",
    ylab = "Population"
  )
  points(years, c_Y1, col = "orange")
  points(years, c_I, col = "darkgreen")
  points(years, c_A, col = "blue")
  points(years, captive_population)
}