#' Function to calculate maximum time weighted average concentrations from
#' kinetic models fitted with mkinfit
#'
#' This function calculates maximum moving window time weighted average
#' concentrations (TWAs) for kinetic models fitted with \code{\link{mkinfit}}.
#' Currently, only calculations for the parent are implemented for the SFO,
#' FOMC, DFOP and HS models, using the analytical formulas given in the PEC
#' soil section of the FOCUS guidance.
#'
#' @aliases max_twa_parent max_twa_sfo max_twa_fomc max_twa_dfop max_twa_hs
#' @param fit An object of class \code{\link{mkinfit}}.
#' @param windows The width of the time windows for which the TWAs should be
#' calculated.
#' @param M0 The initial concentration for which the maximum time weighted
#' average over the decline curve should be calculated. The default is to use
#' a value of 1, which means that a relative maximum time weighted average
#' factor (f_twa) is calculated.
#' @param k The rate constant in the case of SFO kinetics.
#' @param t The width of the time window.
#' @param alpha Parameter of the FOMC model.
#' @param beta Parameter of the FOMC model.
#' @param k1 The first rate constant of the DFOP or the HS kinetics.
#' @param k2 The second rate constant of the DFOP or the HS kinetics.
#' @param g Parameter of the DFOP model.
#' @param tb Parameter of the HS model.
#' @return For \code{max_twa_parent}, a numeric vector, named using the
#' \code{windows} argument. For the other functions, a numeric vector of
#' length one (also known as 'a number').
#' @author Johannes Ranke
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence
#' and Degradation Kinetics from Environmental Fate Studies on Pesticides in
#' EU Registration} Report of the FOCUS Work Group on Degradation Kinetics,
#' EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
#' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
#' @examples
#'
#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
#' max_twa_parent(fit, c(7, 21))
#'
#' @export
max_twa_parent <- function(fit, windows) {
parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
obs_vars <- fit$obs_vars
if (length(obs_vars) > 1) {
warning("Calculation of maximum time weighted average concentrations is",
"currently only implemented for the parent compound using",
"analytical solutions")
}
obs_var <- obs_vars[1]
spec = fit$mkinmod$spec
type = spec[[1]]$type
M0 <- parms.all[paste0(obs_var, "_0")]
if (type == "SFO") {
k_name <- paste0("k_", obs_var)
if (fit$mkinmod$use_of_ff == "min") {
k_name <- paste0(k_name, "_sink")
}
k <- parms.all[k_name]
twafunc <- function(t) {
max_twa_sfo(M0, k, t)
}
}
if (type == "FOMC") {
alpha <- parms.all["alpha"]
beta <- parms.all["beta"]
twafunc <- function(t) {
max_twa_fomc(M0, alpha, beta, t)
}
}
if (type == "DFOP") {
k1 <- parms.all["k1"]
k2 <- parms.all["k2"]
g <- parms.all["g"]
twafunc <- function(t) {
max_twa_dfop(M0, k1, k2, g, t)
}
}
if (type == "HS") {
k1 <- parms.all["k1"]
k2 <- parms.all["k2"]
tb <- parms.all["tb"]
twafunc <- function(t) {
ifelse(t <= tb,
max_twa_sfo(M0, k1, t),
max_twa_hs(M0, k1, k2, tb, t)
)
}
}
if (type %in% c("IORE", "SFORB")) {
stop("Calculation of maximum time weighted average concentrations is currently ",
"not implemented for the ", type, " model.")
}
res <- twafunc(windows)
names(res) <- windows
return(res)
}
#' @rdname max_twa_parent
#' @export
max_twa_sfo <- function(M0 = 1, k, t) {
M0 * (1 - exp(- k * t)) / (k * t)
}
#' @rdname max_twa_parent
#' @export
max_twa_fomc <- function(M0 = 1, alpha, beta, t) {
M0 * (beta)/(t * (1 - alpha)) * ((t/beta + 1)^(1 - alpha) - 1)
}
#' @rdname max_twa_parent
#' @export
max_twa_dfop <- function(M0 = 1, k1, k2, g, t) {
M0/t * ((g/k1) * (1 - exp(- k1 * t)) + ((1 - g)/k2) * (1 - exp(- k2 * t)))
}
#' @rdname max_twa_parent
#' @export
max_twa_hs <- function(M0 = 1, k1, k2, tb, t) {
(M0 / t) * (
(1/k1) * (1 - exp(- k1 * tb)) +
(exp(- k1 * tb) / k2) * (1 - exp(- k2 * (t - tb)))
)
}
