## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4.5,
  dpi       = 120
)

## ----load---------------------------------------------------------------------
library(okBATHTUB)

## ----load_step----------------------------------------------------------------
result <- ok_load(
  inflow_m3yr   = 45e6,    # tributary inflow (m^3/yr)
  tp_inflow_ugl = 120,     # flow-weighted mean TP (ug/L)
  tn_inflow_ugl = 1800     # flow-weighted mean TN (ug/L)
)
print(result)

## ----hydraulics_step----------------------------------------------------------
result <- result |>
  ok_hydraulics(
    surface_area_ha = 890,
    mean_depth_m    = 4.2
  )

cat("Hydraulic residence time:", round(result$data$hydraulic_residence_time_yr, 2), "yr\n")
cat("Areal water load (qs):   ", round(result$data$areal_water_load_myr, 2), "m/yr\n")

## ----retention_step-----------------------------------------------------------
result <- result |> ok_retention()

cat("TP retention coefficient:", round(result$data$tp_retention_coeff, 3),
    "(", result$data$tp_retention_form, ")\n")
cat("TN retention coefficient:", round(result$data$tn_retention_coeff, 3), "\n")

## ----inlake_step--------------------------------------------------------------
result <- result |> ok_inlake()

cat("In-lake TP:    ", round(result$data$tp_inlake_ugl, 1), "ug/L\n")
cat("Chlorophyll-a: ", round(result$data$chla_ugl, 2),    "ug/L\n")
cat("Secchi depth:  ", round(result$data$secchi_m, 2),    "m\n")

## ----tsi_step-----------------------------------------------------------------
result <- result |> ok_tsi()
summary(result)

## ----full_pipeline------------------------------------------------------------
result <- ok_load(
    inflow_m3yr   = 45e6,
    tp_inflow_ugl = 120,
    tn_inflow_ugl = 1800
  ) |>
  ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

summary(result)

## ----reservoir_lookup---------------------------------------------------------
res <- ok_reservoir("Arcadia Lake", exact = TRUE)
res[, c("lake_name", "surface_area_ha", "mean_depth_m", "eco_l3_name",
        "data_quality")]

## ----reservoir_pipeline-------------------------------------------------------
result <- ok_load(
    inflow_m3yr   = 45e6,
    tp_inflow_ugl = 120,
    tn_inflow_ugl = 1800
  ) |>
  ok_hydraulics(
    surface_area_ha = res$surface_area_ha[1],
    mean_depth_m    = res$mean_depth_m[1]
  ) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

cat("Trophic state:", result$data$trophic_state, "\n")
cat("Mean TSI:     ", round(result$data$tsi_mean, 1), "\n")

## ----reservoir_summary--------------------------------------------------------
ok_reservoir_summary()

## ----multi_trib---------------------------------------------------------------
tributaries <- data.frame(
  inflow_m3yr   = c(30e6, 15e6),
  tp_inflow_ugl = c(100,  160),
  tn_inflow_ugl = c(1500, 2400)
)

result_multi <- ok_load_multi(tributaries) |>
  ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
  ok_retention() |>
  ok_inlake()    |>
  ok_tsi()

cat("FW mean TP inflow:", round(result_multi$data$tp_inflow_ugl, 1), "ug/L\n")

## ----tsi_observed-------------------------------------------------------------
obs <- ok_tsi_observed(
  tp_ugl    = 85,
  chla_ugl  = 22,
  secchi_m  = 0.8
)

cat("TSI(TP):    ", round(obs$tsi_tp,     1), "\n")
cat("TSI(Chl-a): ", round(obs$tsi_chla,   1), "\n")
cat("TSI(Secchi):", round(obs$tsi_secchi, 1), "\n")
cat("Mean TSI:   ", round(obs$tsi_mean,   1), "\n")
cat("Class:      ", obs$trophic_state,         "\n")

## ----coeff_comparison---------------------------------------------------------
run_pipeline <- function(coef, eco = NULL) {
  ok_load(
      inflow_m3yr   = 45e6,
      tp_inflow_ugl = 120,
      tn_inflow_ugl = 1800,
      coefficients  = coef,
      ecoregion     = eco
    ) |>
    ok_hydraulics(surface_area_ha = 890, mean_depth_m = 4.2) |>
    ok_retention() |>
    ok_inlake()    |>
    ok_tsi()
}

r_walker  <- run_pipeline("walker")
r_voll    <- run_pipeline("vollenweider")
r_okla_ct <- run_pipeline("oklahoma", "Cross Timbers")

comparison <- data.frame(
  scenario   = c("Walker Model 1",
                 "Vollenweider/Larsen-Mercier",
                 "Oklahoma (Cross Timbers)"),
  tp_inlake  = round(c(r_walker$data$tp_inlake_ugl,
                       r_voll$data$tp_inlake_ugl,
                       r_okla_ct$data$tp_inlake_ugl), 1),
  chla       = round(c(r_walker$data$chla_ugl,
                       r_voll$data$chla_ugl,
                       r_okla_ct$data$chla_ugl), 2),
  tsi_mean   = round(c(r_walker$data$tsi_mean,
                       r_voll$data$tsi_mean,
                       r_okla_ct$data$tsi_mean), 1)
)
print(comparison, row.names = FALSE)