####################
# 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)

simulate_condor_population_inner <- function(rd1, rdI, rdA,
                                             c_Y0, c_Y1, c_I, c_A,
                                             w_Y0, w_Y1, w_I, w_A,
                                             nyears, captive_dc) {
  captive_prod <- if (captive_dc) captive_productivity_rate_dc else captive_productivity_rate
  for (t in (length(c_Y0) + 1):nyears) {
    # 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_prod * (c_A[t] / 2)

    # sum of all population groups
    last_captive_population <- c_Y0[t - 1] + c_Y1[t - 1] + c_I[t - 1] + c_A[t - 1]
    released_to_wild <- if (last_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

    # carrying capacity death rate
    last_wild_population <- w_Y0[t - 1] + w_Y1[t - 1] + w_I[t - 1] + w_A[t - 1]
    rd_carrying <- (last_wild_population / wild_carrying_capacity) # P / M

    w_Y1[t] <- w_Y0[t - 1] - # inflow from year zero
      rd1 * w_Y1[t - 1] # outflow to deaths (death rate)

    w_I[t] <- w_I[t - 1] + # last year's population
      w_Y1[t - 1] + # inflow from year one
      released_to_wild * (1 - rd_carrying) - # inflow from captive
      w_I[t - 1] * wild_immature_to_adult_rate - # outflow to adult
      w_I[t - 1] * rdI # outflow to deaths (death rate)

    w_A[t] <- w_A[t - 1] + # last year's population
      w_I[t - 1] * wild_immature_to_adult_rate - # inflow from immature
      w_A[t - 1] * rdA # outflow to deaths (death rate)

    w_Y0[t] <- wild_productivity_rate * (w_A[t] / 2) *
      (1 - rd_carrying) # deaths due to carrying capacity
  }

  list(
    c_Y0 = c_Y0, c_Y1 = c_Y1, c_I = c_I, c_A = c_A,
    w_Y0 = w_Y0, w_Y1 = w_Y1, w_I = w_I, w_A = w_A
  )
}

plot_populations <- function(c_Y0, c_Y1, c_I, c_A,
                             w_Y0, w_Y1, w_I, w_A) {
  years <- seq(1, simulation_years)

  captive_population <- c_Y0 + c_Y1 + c_I + c_A
  plot(
    years,
    c_Y0,
    ylim = c(0, 150),
    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)

  wild_population <- w_Y0 + w_Y1 + w_I + w_A
  plot(
    years,
    w_Y0,
    ylim = c(0, wild_carrying_capacity),
    col = "red",
    main = "Condor Population in the Wild",
    xlab = "Time (years)",
    ylab = "Population"
  )
  points(years, w_Y1, col = "orange")
  points(years, w_I, col = "darkgreen")
  points(years, w_A, col = "blue")
  points(years, wild_population)

  print(wild_population)
}

# 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, dc = FALSE) {
  # captive population groups
  c_Y0 <- c(30)
  c_Y1 <- c(30)
  c_I <- c(30)
  c_A <- c(30)

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

  list2env(simulate_condor_population_inner(
    rd1, rdI, rdA,
    c_Y0, c_Y1, c_I, c_A,
    w_Y0, w_Y1, w_I, w_A,
    simulation_years, dc
  ), envir = environment())

  plot_populations(
    c_Y0, c_Y1, c_I, c_A,
    w_Y0, w_Y1, w_I, w_A
  )
}

simulate_part_f <- function(rd1, rdI, rdA) {
  # captive population groups
  c_Y0 <- c(30)
  c_Y1 <- c(30)
  c_I <- c(30)
  c_A <- c(30)

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

  sdc <- 5

  # simulate first 5 years
  list2env(simulate_condor_population_inner(
    rd1, rdI, rdA,
    c_Y0, c_Y1, c_I, c_A,
    w_Y0, w_Y1, w_I, w_A,
    sdc, TRUE
  ), envir = environment())

  # release all in captivity
  w_Y0[sdc] <- w_Y0[sdc] + c_Y0[sdc]
  w_Y1[sdc] <- w_Y1[sdc] + c_Y1[sdc]
  w_I[sdc] <- w_I[sdc] + c_I[sdc]
  w_A[sdc] <- w_A[sdc] + c_A[sdc]

  c_Y0[sdc] <- 0
  c_Y1[sdc] <- 0
  c_I[sdc] <- 0
  c_A[sdc] <- 0

  # continue simulation
  list2env(simulate_condor_population_inner(
    rd1, rdI, rdA,
    c_Y0, c_Y1, c_I, c_A,
    w_Y0, w_Y1, w_I, w_A,
    simulation_years, FALSE
  ), envir = environment())

  plot_populations(
    c_Y0, c_Y1, c_I, c_A,
    w_Y0, w_Y1, w_I, w_A
  )
}