213 lines
6.2 KiB
R
213 lines
6.2 KiB
R
####################
|
|
# 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
|
|
)
|
|
}
|