## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4.5,
  message   = FALSE,
  warning   = FALSE
)
library(fb4package)
`%||%` <- function(a, b) if (!is.null(a)) a else b

## ----species------------------------------------------------------------------
data(fish4_parameters)
chinook_key <- grep("tshawytscha", names(fish4_parameters), value = TRUE)[1]
db    <- fish4_parameters[[chinook_key]]
stage <- if ("adult" %in% names(db$life_stages)) "adult" else names(db$life_stages)[1]
sp    <- db$life_stages[[stage]]

cat("Species:", chinook_key, "| Stage:", stage, "\n")
cat("CEQ:", sp$consumption$CEQ,
    " CA:", sp$consumption$CA,
    " CB:", sp$consumption$CB, "\n")

## ----get-p--------------------------------------------------------------------
gp <- function(cat, nm) {
  v <- sp[[cat]][[nm]]
  if (is.null(v)) NA_real_ else as.numeric(v)
}

CEQ <- gp("consumption","CEQ"); CA  <- gp("consumption","CA")
CB  <- gp("consumption","CB");  CQ  <- gp("consumption","CQ")
CTO <- gp("consumption","CTO"); CTM <- gp("consumption","CTM")
CTL <- gp("consumption","CTL"); CK1 <- gp("consumption","CK1")
CK4 <- gp("consumption","CK4")

REQ <- gp("respiration","REQ"); RA  <- gp("respiration","RA")
RB  <- gp("respiration","RB");  RQ  <- gp("respiration","RQ")
RTO <- gp("respiration","RTO"); RTM <- gp("respiration","RTM")
RTL <- gp("respiration","RTL"); RK1 <- gp("respiration","RK1")
RK4 <- gp("respiration","RK4"); RK5 <- gp("respiration","RK5")

ACT  <- gp("activity","ACT");  BACT <- gp("activity","BACT")
SDA  <- gp("sda","SDA")

EGEQ <- gp("egestion","EGEQ"); FA <- gp("egestion","FA")
FB_  <- gp("egestion","FB");   FG <- gp("egestion","FG")
EXEQ <- gp("excretion","EXEQ"); UA <- gp("excretion","UA")
UB   <- gp("excretion","UB");  UG <- gp("excretion","UG")

PREDEDEQ <- gp("predator","PREDEDEQ")

Oxycal <- 13560

## ----temp-coeffs--------------------------------------------------------------
if (!is.na(CEQ) && CEQ == 2) {
  CY <- log(CQ)*(CTM-CTO+2); CZ <- log(CQ)*(CTM-CTO)
  CX <- (CZ^2*(1+(1+40/CY)^0.5)^2)/400
}
if (!is.na(CEQ) && CEQ == 3) {
  CG1 <- (1/(CTO-CQ))*log((0.98*(1-CK1))/(CK1*0.02))
  CG2 <- (1/(CTL-CTM))*log((0.98*(1-CK4))/(CK4*0.02))
}
if (!is.na(REQ) && REQ == 2) {
  RY <- log(RQ)*(RTM-RTO+2); RZ <- log(RQ)*(RTM-RTO)
  RX <- (RZ^2*(1+(1+40/RY)^0.5)^2)/400
}

## ----inputs-------------------------------------------------------------------
set.seed(123)
N_DIAS  <- 180
Init_W  <- 500     # g
p_fixed <- 0.4
ED_ini  <- 4200; ED_fin <- 4600

days     <- 1:N_DIAS
base_T   <- 10 + 4*sin(2*pi*(days-30)/N_DIAS)
Temp_vec <- round(pmax(2, base_T + rnorm(N_DIAS, 0, 0.3)), 2)

prey_prop_mat <- matrix(c(rep(0.7, N_DIAS), rep(0.3, N_DIAS)),
                        nrow = N_DIAS,
                        dimnames = list(NULL, c("Alewife","Inverts")))
prey_ED_mat   <- matrix(c(rep(4900, N_DIAS), rep(2600, N_DIAS)),
                        nrow = N_DIAS,
                        dimnames = list(NULL, c("Alewife","Inverts")))
ind_prey_mat  <- matrix(0, nrow = N_DIAS, ncol = 2,
                        dimnames = list(NULL, c("Alewife","Inverts")))
Pred_E_vec    <- seq(ED_ini, ED_fin, length.out = N_DIAS + 1)

cat(sprintf("%d days | W0 = %.0f g | p = %.2f | T = %.1f–%.1f °C\n",
            N_DIAS, Init_W, p_fixed, min(Temp_vec), max(Temp_vec)))

## ----shiny-fns----------------------------------------------------------------
# --- Temperature functions ---
Cf1T <- function(T) exp(CQ*T)
Cf2T <- function(T) {
  if (T < CTM) { V <- (CTM-T)/(CTM-CTO); ft <- V^CX*exp(CX*(1-V)) } else ft <- 0
  max(ft, 0)
}
Cf3T <- function(T) {
  L1 <- exp(CG1*(T-CQ)); KA <- (CK1*L1)/(1+CK1*(L1-1))
  L2 <- exp(CG2*(CTL-T)); KB <- (CK4*L2)/(1+CK4*(L2-1))
  KA*KB
}
Cf4T <- function(T) max(exp(CQ*T + CK1*T^2 + CK4*T^3), 0)

Rf1T    <- function(T) exp(RQ*T)
RACTf1T <- function(W,T) {
  VEL <- if (T <= RTL) ACT*W^RK4*exp(BACT*T) else RK1*W^RK4*exp(RK5*T)
  exp(RTO*VEL)
}
Rf2T <- function(T) {
  if (T < RTM) { V <- (RTM-T)/(RTM-RTO); ft <- V^RX*exp(RX*(1-V)) } else ft <- 1e-6
  max(ft, 1e-6)
}

# --- Bioenergetic functions ---
consumo_fn <- function(T, W, p) {
  Cmax <- CA*W^CB
  ft   <- switch(as.character(CEQ),
                 "1" = Cf1T(T), "2" = Cf2T(T),
                 "3" = Cf3T(T), "4" = Cf4T(T))
  Cmax*p*ft
}

respiracion_fn <- function(T, W) {
  Rmax <- RA*W^RB
  if (REQ == 1) { ft <- Rf1T(T); ACT_v <- RACTf1T(W,T) } else { ft <- Rf2T(T); ACT_v <- ACT }
  Rmax*ft*ACT_v
}

egestion_fn <- function(C, T, p, i) {
  switch(as.character(EGEQ),
    "1" = FA*C,
    "2" = FA*(T^FB_)*exp(FG*p)*C,
    "3" = { PE  <- FA*(T^FB_)*exp(FG*p)
            PFF <- sum(ind_prey_mat[i,]*prey_prop_mat[i,])
            PF  <- ((PE-0.1)/0.9)*(1-PFF)+PFF; PF*C },
    "4" = FA*(T^FB_)*C)
}

excrecion_fn <- function(C, Eg, T, p) {
  switch(as.character(EXEQ),
    "1" = UA*(C-Eg),
    "2" =, "3" = UA*(T^UB)*exp(UG*p)*(C-Eg),
    "4" = UA*(T^UB)*(C-Eg))
}

# --- Daily growth loop ---
grow_shiny_daily <- function(Temperature, W, mode = "p_value", mode_value = 0.4,
                             prey_pm, prey_em, Fin) {
  out <- vector("list", Fin)
  for (i in seq_len(Fin)) {
    Pred_E_i  <- Pred_E_vec[min(i,   length(Pred_E_vec))]
    Pred_E_ip1 <- Pred_E_vec[min(i+1, length(Pred_E_vec))]
    mean_prey <- sum(prey_pm[i,] * prey_em[i,])

    if (mode == "p_value") {
      p       <- mode_value
      Cons_gg <- consumo_fn(Temperature[i], W, p)
    } else if (mode == "Ration") {
      Cons_gg <- mode_value
      Cmax    <- consumo_fn(Temperature[i], W, 1.0)
      p       <- if (Cmax > 0) Cons_gg/Cmax else 0
    } else if (mode == "Ration_prey") {
      Cons_gg <- mode_value / W
      Cmax    <- consumo_fn(Temperature[i], W, 1.0)
      p       <- if (Cmax > 0) Cons_gg/Cmax else 0
    }

    Cons <- Cons_gg * mean_prey
    Eg   <- egestion_fn(Cons, Temperature[i], p, i)
    Ex   <- excrecion_fn(Cons, Eg, Temperature[i], p)
    S    <- SDA*(Cons-Eg)
    Res  <- respiracion_fn(Temperature[i], W)*Oxycal

    G      <- Cons-(Res+Eg+Ex+S)
    egain  <- G*W
    finalwt <- (egain + Pred_E_i*W) / Pred_E_ip1
    if (is.na(finalwt) || finalwt < 0) finalwt <- 0

    out[[i]] <- list(
      Day                  = i,
      Starting_Weight      = W,
      Weight               = finalwt,
      Weight_change        = finalwt - W,
      Temperature          = Temperature[i],
      Mean_Prey_Energy_J_g = mean_prey,
      Consumption_gg       = Cons_gg,
      Consumption_energy   = Cons,
      Respiration          = Res,
      Egestion             = Eg,
      Excretion            = Ex,
      SDA                  = S,
      Net_energy           = G,
      Energy_density       = Pred_E_ip1,
      Cons_Alewife_g       = Cons_gg*W*prey_pm[i,"Alewife"],
      Cons_Inverts_g       = Cons_gg*W*prey_pm[i,"Inverts"]
    )
    W <- finalwt
  }
  as.data.frame(do.call(rbind, lapply(out, as.data.frame)))
}

# --- Binary search (Weight or Consumption target) ---
fit_p_shiny <- function(IW, target, fit_type, Temperature,
                        prey_pm, prey_em, Fin,
                        W_tol = 0.001, max_iter = 100) {
  p_min <- 0; p_max <- 5; p <- 0.5

  sim_val <- function(p_try) {
    df <- grow_shiny_daily(Temperature, IW, "p_value", p_try, prey_pm, prey_em, Fin)
    if (fit_type == "Weight") tail(df$Weight, 1)
    else sum(df$Consumption_gg * df$Starting_Weight)
  }

  val <- sim_val(p)
  for (iter in seq_len(max_iter)) {
    if (abs(val - target) <= W_tol) break
    if (val > target) p_max <- p else p_min <- p
    p <- (p_min + p_max)/2
    val <- sim_val(p)
  }
  list(p = p, converged = (abs(val-target) <= W_tol), final_val = val)
}

## ----pkg-object---------------------------------------------------------------
sp_pkg <- sp
sp_pkg$predator$PREDEDEQ <- 1L
sp_pkg$predator$ED_ini   <- ED_ini
sp_pkg$predator$ED_end   <- ED_fin

temp_df   <- data.frame(Day = days, Temperature = Temp_vec)
diet_df   <- data.frame(Day = days, Alewife = prey_prop_mat[,"Alewife"],
                                     Inverts = prey_prop_mat[,"Inverts"])
energy_df <- data.frame(Day = days, Alewife = prey_ED_mat[,"Alewife"],
                                     Inverts = prey_ED_mat[,"Inverts"])

bio <- Bioenergetic(
  species_params      = sp_pkg,
  species_info        = list(scientific_name = "Oncorhynchus tshawytscha",
                              common_name = "Chinook", life_stage = "adult"),
  environmental_data  = list(temperature = temp_df),
  diet_data           = list(proportions = diet_df,
                             prey_names  = c("Alewife","Inverts"),
                             energies    = energy_df),
  simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
)

## ----compare-fn---------------------------------------------------------------
compare_outputs <- function(shiny_df, pkg_df, method_label,
  cols = c("Starting_Weight","Weight","Weight_change",
           "Mean_Prey_Energy_J_g","Consumption_gg","Consumption_energy",
           "Respiration","Egestion","Excretion","SDA",
           "Net_energy","Energy_density","Cons_Alewife_g","Cons_Inverts_g")) {

  cols_ok <- intersect(cols, intersect(names(shiny_df), names(pkg_df)))

  res <- lapply(cols_ok, function(col) {
    s <- as.numeric(shiny_df[[col]])
    p <- as.numeric(pkg_df[[col]])
    if (all(is.na(s)) || all(is.na(p))) return(NULL)
    d    <- p - s
    rmse <- sqrt(mean(d^2, na.rm = TRUE))
    r2   <- if (var(s, na.rm=TRUE) > 0) cor(s, p, use="complete.obs")^2 else NA
    data.frame(method = method_label, col = col,
               max_abs_diff = max(abs(d), na.rm=TRUE),
               rmse = rmse, r2 = r2)
  })
  do.call(rbind, Filter(Negate(is.null), res))
}

## ----m1, cache=TRUE-----------------------------------------------------------
shiny1 <- grow_shiny_daily(Temp_vec, Init_W, "p_value", p_fixed,
                           prey_prop_mat, prey_ED_mat, N_DIAS)
pkg1   <- run_fb4(bio, fit_to = "p_value", fit_value = p_fixed,
                  strategy = "direct", verbose = FALSE)$daily_output

m1 <- compare_outputs(shiny1, pkg1, "1 — p-value (direct)")

## ----m2, cache=TRUE-----------------------------------------------------------
target_W <- tail(shiny1$Weight, 1)
fit2     <- fit_p_shiny(Init_W, target_W, "Weight",
                        Temp_vec, prey_prop_mat, prey_ED_mat, N_DIAS)
shiny2   <- grow_shiny_daily(Temp_vec, Init_W, "p_value", fit2$p,
                             prey_prop_mat, prey_ED_mat, N_DIAS)
pkg2     <- run_fb4(bio, fit_to = "Weight", fit_value = target_W,
                    strategy = "binary_search", verbose = FALSE)$daily_output

m2 <- compare_outputs(shiny2, pkg2, "2 — Weight (binary search)")

## ----m3, cache=TRUE-----------------------------------------------------------
target_C <- sum(shiny1$Consumption_gg * shiny1$Starting_Weight)
fit3     <- fit_p_shiny(Init_W, target_C, "Consumption",
                        Temp_vec, prey_prop_mat, prey_ED_mat, N_DIAS)
shiny3   <- grow_shiny_daily(Temp_vec, Init_W, "p_value", fit3$p,
                             prey_prop_mat, prey_ED_mat, N_DIAS)
pkg3     <- run_fb4(bio, fit_to = "Consumption", fit_value = target_C,
                    strategy = "binary_search", verbose = FALSE)$daily_output

m3 <- compare_outputs(shiny3, pkg3, "3 — Consumption (binary search)")

## ----m4, cache=TRUE-----------------------------------------------------------
Ration_pct  <- 1.5       # % body weight per day
Ration_frac <- Ration_pct / 100
shiny4 <- grow_shiny_daily(Temp_vec, Init_W, "Ration", Ration_frac,
                           prey_prop_mat, prey_ED_mat, N_DIAS)
pkg4   <- run_fb4(bio, fit_to = "Ration", fit_value = Ration_pct,
                  strategy = "direct_ration_percent", verbose = FALSE)$daily_output

m4 <- compare_outputs(shiny4, pkg4, "4 — Ration % body weight (direct)")

## ----m5, cache=TRUE-----------------------------------------------------------
Ration_g <- 8.0     # g/day
shiny5   <- grow_shiny_daily(Temp_vec, Init_W, "Ration_prey", Ration_g,
                             prey_prop_mat, prey_ED_mat, N_DIAS)
pkg5     <- run_fb4(bio, fit_to = "Ration_prey", fit_value = Ration_g,
                    strategy = "direct", verbose = FALSE)$daily_output

m5 <- compare_outputs(shiny5, pkg5, "5 — Ration g/day (direct)")

## ----summary-table------------------------------------------------------------
all_metrics <- rbind(m1, m2, m3, m4, m5)

# Worst RMSE per method (excluding Energy_density — known reporting offset)
core_cols <- setdiff(unique(all_metrics$col), "Energy_density")
worst <- do.call(rbind, lapply(split(all_metrics[all_metrics$col %in% core_cols, ],
                                     all_metrics[all_metrics$col %in% core_cols, "method"]),
                               function(x) x[which.max(x$rmse), ]))

knitr::kable(
  worst[, c("method", "col", "max_abs_diff", "rmse", "r2")],
  digits   = c(0, 0, 2e-10, 2e-10, 8),
  format.args = list(scientific = TRUE),
  col.names   = c("Method", "Worst column", "Max |diff|", "RMSE", "R²"),
  caption     = "Worst-case RMSE per method (core bioenergetic variables)"
)

## ----plot-weight, fig.cap="Daily weight: fb4package (solid) vs FB4-Shiny (dashed). Lines overlap completely."----
plot(shiny1$Day, shiny1$Weight, type = "l", lty = 2, col = "steelblue",
     xlab = "Day", ylab = "Weight (g)", main = "Method 1 — p-value direct")
lines(pkg1$Day, pkg1$Weight, lty = 1, col = "firebrick", lwd = 2)
legend("topleft", legend = c("fb4package","FB4-Shiny"),
       col = c("firebrick","steelblue"), lty = c(1,2), lwd = c(2,1), bty = "n")

## ----plot-resid, fig.cap="Daily weight difference (package − Shiny). Scale is sub-femtogram for Method 1."----
par(mfrow = c(1, 2))

diff1 <- pkg1$Weight - shiny1$Weight
plot(diff1, type = "l", col = "firebrick",
     main = "Method 1 residuals", xlab = "Day", ylab = "Δ Weight (g)")
abline(h = 0, lty = 2)

diff2 <- pkg2$Weight - shiny2$Weight
plot(diff2, type = "l", col = "steelblue",
     main = "Method 2 residuals", xlab = "Day", ylab = "Δ Weight (g)")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))

## ----nutrient-val, cache=TRUE-------------------------------------------------
prey_n <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.025,0.025), Inverts = c(0.018,0.018))
prey_p <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.004,0.004), Inverts = c(0.003,0.003))
pred_n <- data.frame(Day = c(1,N_DIAS), value = c(0.030,0.030))
pred_p <- data.frame(Day = c(1,N_DIAS), value = c(0.004,0.004))
n_ae   <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.80,0.80), Inverts = c(0.75,0.75))
p_ae   <- data.frame(Day = c(1,N_DIAS), Alewife = c(0.60,0.60), Inverts = c(0.55,0.55))

bio_nut <- Bioenergetic(
  species_params      = sp_pkg,
  species_info        = list(scientific_name = "Oncorhynchus tshawytscha"),
  environmental_data  = list(temperature = temp_df),
  diet_data           = list(proportions = diet_df, prey_names = c("Alewife","Inverts"),
                             energies = energy_df),
  nutrient_data       = list(prey_n_concentrations    = prey_n,
                             prey_p_concentrations    = prey_p,
                             predator_n_concentration = pred_n,
                             predator_p_concentration = pred_p,
                             n_assimilation_efficiency = n_ae,
                             p_assimilation_efficiency = p_ae),
  simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
)

nut <- run_fb4(bio_nut, fit_to = "p_value", fit_value = p_fixed,
               strategy = "direct", verbose = FALSE)$daily_output

N_err <- with(nut, abs(Nitrogen_consumed_g -
                       (Nitrogen_growth_g + Nitrogen_excretion_g + Nitrogen_egestion_g)))
P_err <- with(nut, abs(Phosphorus_consumed_g -
                       (Phosphorus_growth_g + Phosphorus_excretion_g + Phosphorus_egestion_g)))

data.frame(
  Check   = c("N mass balance (max error, g)", "P mass balance (max error, g)",
              "N consumed ≥ 0", "P consumed ≥ 0"),
  Result  = c(format(max(N_err, na.rm=TRUE), scientific=TRUE),
              format(max(P_err, na.rm=TRUE), scientific=TRUE),
              all(nut$Nitrogen_consumed_g   >= 0, na.rm=TRUE),
              all(nut$Phosphorus_consumed_g >= 0, na.rm=TRUE))
) |> knitr::kable(col.names = c("Check", "Result"),
                  caption = "Nutrient mass balance checks")

## ----plot-nutrients, fig.cap="Daily nitrogen and phosphorus fluxes over 180 days."----
par(mfrow = c(1, 2))
plot(nut$Day, nut$Nitrogen_consumed_g, type = "l", col = "steelblue",
     xlab = "Day", ylab = "g N / day", main = "Nitrogen fluxes")
lines(nut$Day, nut$Nitrogen_growth_g,    col = "darkgreen")
lines(nut$Day, nut$Nitrogen_excretion_g, col = "firebrick", lty = 2)
lines(nut$Day, nut$Nitrogen_egestion_g,  col = "orange",    lty = 2)
legend("topright", legend = c("Consumed","Growth","Excretion","Egestion"),
       col = c("steelblue","darkgreen","firebrick","orange"),
       lty = c(1,1,2,2), bty = "n", cex = 0.8)

plot(nut$Day, nut$Phosphorus_consumed_g, type = "l", col = "steelblue",
     xlab = "Day", ylab = "g P / day", main = "Phosphorus fluxes")
lines(nut$Day, nut$Phosphorus_growth_g,    col = "darkgreen")
lines(nut$Day, nut$Phosphorus_excretion_g, col = "firebrick", lty = 2)
lines(nut$Day, nut$Phosphorus_egestion_g,  col = "orange",    lty = 2)
par(mfrow = c(1, 1))

## ----contam-val, cache=TRUE---------------------------------------------------
cont_prey <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.05,0.05), Inverts=c(0.02,0.02))
cont_ae   <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.80,0.80), Inverts=c(0.70,0.70))
cont_te   <- data.frame(Day=c(1,N_DIAS), Alewife=c(0.85,0.85), Inverts=c(0.75,0.75))

run_conteq <- function(CONTEQ_num, extra = list()) {
  cd <- c(list(CONTEQ = CONTEQ_num, initial_concentration = 0.10,
               prey_concentrations = cont_prey,
               transfer_efficiency = cont_te,
               assimilation_efficiency = cont_ae), extra)
  b <- Bioenergetic(
    species_params      = sp_pkg,
    species_info        = list(scientific_name = "Oncorhynchus tshawytscha"),
    environmental_data  = list(temperature = temp_df),
    diet_data           = list(proportions = diet_df, prey_names = c("Alewife","Inverts"),
                               energies = energy_df),
    contaminant_data    = cd,
    simulation_settings = list(initial_weight = Init_W, duration = N_DIAS)
  )
  run_fb4(b, fit_to = "p_value", fit_value = p_fixed,
          strategy = "direct", verbose = FALSE)$daily_output
}

df_c1 <- run_conteq(1)
df_c2 <- run_conteq(2)
df_c3 <- run_conteq(3, list(gill_efficiency = 0.5, fish_water_partition = 5000,
                             water_concentration = 0.0001, dissolved_fraction = 0.9,
                             do_saturation = 0.9))

cont_checks <- data.frame(
  Check = c(
    "CONTEQ 1 — burden ≥ 0",
    "CONTEQ 1 — clearance = 0 (pure accumulation)",
    "CONTEQ 1 — burden monotone",
    "CONTEQ 2 — burden ≥ 0",
    "CONTEQ 2 — elimination active (>50% of days)",
    "CONTEQ 3 — burden ≥ 0",
    "CONTEQ 3 — elimination active (>50% of days)"
  ),
  Pass = c(
    all(df_c1$Contaminant_burden_ug >= 0, na.rm=TRUE),
    all(df_c1$Contaminant_clearance_ug_d == 0, na.rm=TRUE),
    all(diff(df_c1$Contaminant_burden_ug) >= -1e-10, na.rm=TRUE),
    all(df_c2$Contaminant_burden_ug >= 0, na.rm=TRUE),
    mean(df_c2$Contaminant_clearance_ug_d > 0, na.rm=TRUE) > 0.5,
    all(df_c3$Contaminant_burden_ug >= 0, na.rm=TRUE),
    mean(df_c3$Contaminant_clearance_ug_d > 0, na.rm=TRUE) > 0.5
  )
)

knitr::kable(cont_checks, caption = "Contaminant sub-model checks")

## ----plot-contam, fig.cap="Contaminant body burden over 180 days for three CONTEQ equations."----
plot(df_c1$Day, df_c1$Contaminant_burden_ug, type = "l",
     col = "firebrick", lwd = 2,
     xlab = "Day", ylab = "Burden (μg)", main = "Contaminant body burden by CONTEQ")
lines(df_c2$Day, df_c2$Contaminant_burden_ug, col = "steelblue", lwd = 2)
lines(df_c3$Day, df_c3$Contaminant_burden_ug, col = "darkgreen", lwd = 2)
legend("topleft", legend = c("CONTEQ 1 (no clearance)",
                             "CONTEQ 2 (T/W-dependent)",
                             "CONTEQ 3 (Arnot & Gobas)"),
       col = c("firebrick","steelblue","darkgreen"), lwd = 2, bty = "n")