Support units in UK drainage calculations

Also, address some R CMD check notes

4 files changed, 37 insertions, 19 deletions

diff --git a/R/PEC_sw_drainage_UK.R b/R/PEC_sw_drainage_UK.R
index 4b3111e..599e15e 100644
--- a/R/PEC_sw_drainage_UK.R
+++ b/R/PEC_sw_drainage_UK.R
@@ -4,11 +4,14 @@
 #' published on the CRC website. Degradation before the start of the drainage period is taken into account if
 #' `latest_application` is specified and the degradation parameters are given either as a `soil_DT50` or a `model`.
 #'
-#' @param rate Application rate in g/ha
+#' @param rate Application rate in g/ha or with a compatible unit specified
+#' with the units package
 #' @param interception The fraction of the application rate that does not reach the soil
-#' @param Koc The sorption coefficient normalised to organic carbon in L/kg
+#' @param Koc The sorption coefficient normalised to organic carbon in L/kg or a unit specified
+#' with the units package
 #' @param latest_application Latest application date, formatted as e.g. "01 July"
-#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used
+#' @param soil_DT50 Soil degradation half-life, if SFO kinetics are to be used, in
+#' days or a time unit specified with the units package
 #' @param model The soil degradation model to be used. Either one of "FOMC",
 #'   "DFOP", "HS", or "IORE", or an mkinmod object
 #' @param model_parms A named numeric vector containing the model parameters
@@ -27,27 +30,41 @@
 #' PEC_sw_drainage_UK(150, Koc = 100)
 #' PEC_sw_drainage_UK(60, interception = 0.5, Koc = 550,
 #'   latest_application = "01 July", soil_DT50 = 200)
-
-PEC_sw_drainage_UK <- function(rate, interception = 0, Koc,
-                                   latest_application = NULL, soil_DT50 = NULL,
-                                   model = NULL, model_parms = NULL)
+PEC_sw_drainage_UK <- function(rate, 
+  interception = 0, Koc,
+  latest_application = NULL, soil_DT50 = NULL,
+  model = NULL, model_parms = NULL)
 {
-  percentage_lost <- SSLRC_mobility_classification(Koc)[[2]]
-  amount_available <- rate * (1 - interception) # g/ha
+  # Set default units if not specified and convert to units used in the calculations
+  if (!inherits(rate, "units")) rate <- set_units(rate, "g/ha")
+  rate_g_ha <- as.numeric(set_units(rate, "g/ha"))
+
+  if (!inherits(Koc, "units")) Koc <- set_units(Koc, "L/kg")
+  Koc_L_kg <- as.numeric(set_units(Koc, "L/kg"))
+
+  if (!missing(soil_DT50)) {
+    if (!inherits(soil_DT50, "units")) {
+      soil_DT50_d <- soil_DT50
+    }
+    soil_DT50_d <- as.numeric(set_units(soil_DT50, "d"))
+  }
+
+  percentage_lost <- SSLRC_mobility_classification(Koc_L_kg)[[2]]
+  amount_available <- rate_g_ha * (1 - interception) # amount in g for 1 ha
 
   if (!missing(latest_application)) {
     lct <- Sys.getlocale("LC_TIME")
     tmp <- Sys.setlocale("LC_TIME", "C")
-    if (latest_application == "29 February") {
+    if (latest_application == "29 February") { # Use a leap year
       ref_year <- 2000
-    } else { ref_year <- 1999}
+    } else { ref_year <- 1999} # Use a non-leap year
     latest <- as.Date(paste(latest_application, ref_year), "%d %b %Y")
     if (is.na(latest)) stop("Please specify the latest application in the format '%d %b', e.g. '01 July'")
     tmp <- Sys.setlocale("LC_TIME", lct)
     degradation_time <- as.numeric(difftime(as.Date(paste0(ref_year,"-10-01")), units = "days", latest))
     if (degradation_time > 0) {
       if (!missing(soil_DT50)) {
-        k = log(2)/soil_DT50
+        k = log(2)/soil_DT50_d
         as.Date(paste(latest_application, "1999"), "%d %B %Y")
 
         amount_available <- amount_available * exp(-k * degradation_time)
@@ -62,6 +79,6 @@ PEC_sw_drainage_UK <- function(rate, interception = 0, Koc,
   } 
 
   volume = 130000 # L/ha
-  PEC = 1e6 * (percentage_lost/100) * amount_available / volume
+  PEC = set_units(1e6 * (percentage_lost/100) * amount_available / volume, "\u00B5g/L")
   return(PEC)
 }
diff --git a/R/PEC_sw_exposit.R b/R/PEC_sw_exposit.R
index 282a1ac..a157189 100644
--- a/R/PEC_sw_exposit.R
+++ b/R/PEC_sw_exposit.R
@@ -52,7 +52,8 @@
 #'
 #' @importFrom units as_units set_units drop_units
 #' @importFrom dplyr across mutate
-#' @param rate The application rate in g/ha
+#' @param rate Application rate in g/ha or with a compatible unit specified
+#' with the units package
 #' @param interception The fraction intercepted by the crop
 #' @param Koc The sorption coefficient to soil organic carbon
 #' @param DT50 The soil half-life in days
diff --git a/R/TOXSWA_cwa.R b/R/TOXSWA_cwa.R
index 310029e..b2f7619 100644
--- a/R/TOXSWA_cwa.R
+++ b/R/TOXSWA_cwa.R
@@ -147,10 +147,10 @@ plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t
 #' and some associated statistics.  like maximum moving window average
 #' concentrations, and dataframes holding the events exceeding specified
 #' thresholds.  Usually, an instance of this class will be generated
-#' by \code{\link{read.TOXSWA_cwa}}.
+#' by [read.TOXSWA_cwa].
 #'
 #' @export
-#' @format An \code{\link{R6Class}} generator object.
+#' @format An [R6::R6Class] generator object.
 #' @field filename Length one character vector holding the filename.
 #' @field basedir Length one character vector holding the directory where the file came from.
 #' @field zipfile If not null, giving the path to the zip file from which the file was read.
diff --git a/R/twa.R b/R/twa.R
index 2018dfa..c756f58 100644
--- a/R/twa.R
+++ b/R/twa.R
@@ -3,9 +3,9 @@
 #' @param x When numeric, this is the half-life to be used for an exponential
 #'   decline. When a character string specifying a parent decline model is given
 #'   e.g. \code{FOMC}, \code{parms} must contain the corresponding parameters.
-#'   If x is an \code{\link{mkinfit}} object, the decline is calculated from this
+#'   If x is an [mkin::mkinfit] object, the decline is calculated from this
 #'   object.
-#' @param ini The initial amount. If x is an \code{\link{mkinfit}} object, and
+#' @param ini The initial amount. If x is an [mkin::mkinfit] object, and
 #'   ini is 'model', the fitted initial concentrations are used. Otherwise, ini
 #'   must be numeric. If it has length one, it is used for the parent and
 #'   initial values of metabolites are zero, otherwise, it must give values for
@@ -246,7 +246,7 @@ twa.one_box <- function(x, window = 21)
 #' \code{\link{plot.one_box}} using the window size for the argument
 #' \code{max_twa}.
 #'
-#' The method working directly on fitted \code{\link{mkinfit}} objects uses the
+#' The method working directly on fitted [mkin::mkinfit]] objects uses the
 #' equations given in the PEC soil section of the FOCUS guidance and is restricted
 #' SFO, FOMC and DFOP models and to the parent compound
 #' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and