## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup,echo = FALSE-------------------------------------------------------
library(DMSTAr)

## ----Fc calc, echo = FALSE,fig.align = "center",fig.cap = "Graphical representation of $F_c$.",fig.width = 6.5,fig.height = 4----
## DMSTA uses Michaelis- Menten / half saturation


Fc <- function(C, Chalf) {
  C / (C + Chalf)
}


C  <- seq(0, 200, by = 10)  # TP concentration (µg/L)
Chalf <- 50                 # typical DMSTA value
Fc_dat <- data.frame(C = C,
                     Fc = Fc(C, Chalf))

oldpar <- par(no.readonly = TRUE)
par(mar=c(2,3,1,1.5),oma = c(2,2,0.5,0.5))
ylim.val <- c(0,1);by.y <- 0.2;ymaj <- seq(ylim.val[1],ylim.val[2],by.y)
xlim.val <- c(0,200);by.x <- 50;xmaj <- seq(xlim.val[1],xlim.val[2],by.x)

plot(Fc~C, Fc_dat,type = "n",ylim=ylim.val,xlim=xlim.val,las=1,ann=F)
abline(h=ymaj,v=xmaj,lty=3,col="grey",lwd=0.5)
lines(Fc~C, Fc_dat,lwd=2)
abline(h = c(0, 1), lty = 2,lwd=2,col="red")
mtext(side=3,adj=0.01,expression("Concentration-dependent scale factor " * F[c]))
mtext(side=1,line=2.5,expression("Water-column TP conc. ("*mu*"g"~"L"^"-1"*")"))
mtext(side=2,line=2.75,expression(F[c]),cex=1.25)

par(oldpar)

## ----Fz calc, echo = FALSE,fig.align = "center",fig.cap = "Graphical representation of $F_z$.",fig.width = 6.5,fig.height = 4----
Fz <- function(z, Z1, Z2, Z3) {
  # Based on DMSTA2 calibration
  ifelse(
    z < Z1, z / Z1,
    ifelse(
      z <= Z2, 1,
      ifelse(
        z < Z3, (Z3 - z) / (Z3 - Z2),
        0
      )
    )
  )
}

z  <- seq(0, 250, by = 1)  # water depth (cm)

Z1 <- 40    # saturated uptake depth
Z2 <- 100   # lower penalty depth
Z3 <- 200   # upper penalty depth

Fz_dat <- data.frame(Z = z,
                     Fz = Fz(z, Z1, Z2, Z3))
oldpar <- par(no.readonly = TRUE)
par(mar=c(2,3,1,1.5),oma = c(2,2,0.5,0.5))
ylim.val <- c(0,1);by.y <- 0.2;ymaj <- seq(ylim.val[1],ylim.val[2],by.y)
xlim.val <- c(0,200);by.x <- 50;xmaj <- seq(xlim.val[1],xlim.val[2],by.x)

plot(Fz~Z, Fz_dat,type = "n",ylim=ylim.val,xlim=xlim.val,las=1,ann=F)
abline(h=ymaj,v=xmaj,lty=3,col="grey",lwd=0.5)
lines(Fz~Z, Fz_dat, lwd=2)
abline(h = c(0, 1), lty = 2,lwd=2,col="red")
mtext(side=3,adj=0.01,expression("Depth-dependent scale factor " * F[z]))
mtext(side=1,line=2.5,"Water Depth (cm)")
mtext(side=2,line=2.75,expression(F[z]),cex=1.25)
par(oldpar)
# Fz_dat$Fz_EAV <- with(Fz_dat,pmin(1,Z/40)); #DMSTA1 
# lines(Fz_EAV~Z, Fz_dat,lty=2 ,lwd=2,col="blue")

## ----Pmodel1------------------------------------------------------------------
DMSTAr:::compute_DMSTA_kvals(C1000 = 22, Cstar = 3, Ks = 16)


## ----PModel2------------------------------------------------------------------
DMSTAr:::compute_DMSTA_kvals(C1000 = 22, Cstar = -3, Ks = 16)


## ----kinetics-----------------------------------------------------------------
pparams <- list(
  # shared / STA
  C1000 = 22, Cstar = 3, Ks_per_yr = 16,
  Z1 = 40, Z2 = 100, Z3 = 200,
  K2Coef1 = 0.1, Chalf = 50, SeasonalFactor = 1,
  # PSTA
  Ytrans = 1, Ysigma = 1, C1000_2 = 50, ks_2 = 20, zh_2 = 10,
  # RES
  k_depth_penalty = 0.5,
  DutyCycle = 0.95
)

kin_out <- build_P_kin_slots(
  mods = c("STA", "PSTA", "RES"),
  pparams = pparams,
  Dpy = 365.25
)
kin_out