Finish the Step 1 calculator including tests

Some cleaning up. PELMO facilities do not currently work at my end,
as I have no working wine installation on this computer

4 files changed, 113 insertions, 27 deletions

diff --git a/R/FOCUS_GW_scenarios_2012.R b/R/FOCUS_GW_scenarios_2012.R
index e26ea50..33c5293 100644
--- a/R/FOCUS_GW_scenarios_2012.R
+++ b/R/FOCUS_GW_scenarios_2012.R
@@ -1,4 +1,4 @@
-#' A very small subset of the FOCUS Groundwater scenario defitions
+#' A very small subset of the FOCUS Groundwater scenario definitions
 #' 
 #' Currently, only scenario names with acronyms and a small subset of the soil definitions are provided. The 
 #' soil definitions are from page 46ff. from FOCUS (2012).
diff --git a/R/FOCUS_Step_12_scenarios.R b/R/FOCUS_Step_12_scenarios.R
index cb0654a..f28f03d 100644
--- a/R/FOCUS_Step_12_scenarios.R
+++ b/R/FOCUS_Step_12_scenarios.R
@@ -1,3 +1,6 @@
+# Register global variables
+if(getRversion() >= '2.15.1') utils::globalVariables("FOCUS_Step_12_scenarios")
+
 #' Step 1/2 scenario data as distributed with the FOCUS Step 1/2 calculator
 #'
 #' The data were extracted from the scenario.txt file using the R code shown below.
diff --git a/R/PEC_sw_focus.R b/R/PEC_sw_focus.R
index 5f3b3d2..21ed514 100644
--- a/R/PEC_sw_focus.R
+++ b/R/PEC_sw_focus.R
@@ -1,37 +1,75 @@
 #' Calculate FOCUS Step 1 PEC surface water
 #'
-#' This is an attempt to reimplement the FOCUS Step 1 and 2 calculator authored
-#' by Michael Klein. The Step 1 and 2 scenario assumes an area ratio of 10:1
-#' between field and waterbody, and a water depth of 30 cm.
-#' I did not (yet) implement the TWA formulas for times later than day 1, as I
-#' did not understand them right away.
-#' Also, Step 2 is not implemented (yet).
+#' This is reimplementation of Step 1 of the FOCUS Step 1 and 2 calculator
+#' version 3.2, authored by Michael Klein. Note that results for multiple
+#' applications should be compared to the corresponding results for a
+#' single application. At current, this is not done automatically in
+#' this implementation.
 #'
+#' @importFrom utils read.table
+#' @references FOCUS (2014) Generic guidance for Surface Water Scenarios (version 1.4).
+#'  FOrum for the Co-ordination of pesticde fate models and their USe.
+#'  http://esdac.jrc.ec.europa.eu/public_path/projects_data/focus/sw/docs/Generic%20FOCUS_SWS_vc1.4.pdf
+#' @references Website of the Steps 1 and 2 calculator at the Joint Research
+#'   Center of the European Union:
+#'   http://esdac.jrc.ec.europa.eu/projects/stepsonetwo
+#' @note The formulas for input to the waterbody via runoff/drainage of the
+#'   parent and subsequent formation of the metabolite in water is not
+#'   documented in the model description coming with the calculator
+#' @note Step 2 is not implemented
 #' @export
 #' @param parent A list containing substance specific parameters
 #' @param rate The application rate in g/ha. Overriden when
 #'   applications are given explicitly
 #' @param n The number of applications
 #' @param i The application interval
-#' @param applications A dataframe containing times and amounts of each application
-#' @param step At the moment, only Step 1 is implemented
+#' @param met A list containing metabolite specific parameters. If not NULL, 
+#'   the PEC is calculated for this compound, not the parent.
+#' @param f_drift The fraction of the application rate reaching the waterbody
+#'   via drift. If NA, this is derived from the scenario name and the number
+#'   of applications via the drift data defined by the 
+#'   \code{\link{FOCUS_Step_12_scenarios}}
+#' @param f_rd The fraction of the amount applied reaching the waterbody via
+#'   runoff/drainage. At Step 1, it is assumed to be 10%, be it the
+#'   parent or a metabolite
+#' @param scenario The name of the scenario. Must be one of the scenario
+#'   names given in \code{\link{FOCUS_Step_12_scenarios}}
 #' @examples
+#' # Parent only
 #' dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
 #' PEC_sw_focus(dummy_1, 3000, f_drift = 0)
+#'
+#' # Metabolite
+#' new_dummy <- chent_focus_sw(mw = 250, Koc = 100)
+#' M1 <- chent_focus_sw(mw = 100, cwsat = 100, DT50_ws = 100, Koc = 50, max_ws = 0, max_soil = 0.5)
+#' PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter", met = M1)
 PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
-  applications = data.frame(time = seq(0, 0 + n * i, length.out = n),
-                            amount = rate),
   met = NULL,
-  step = 1,
-  f_drift = 0.02759, f_rd = 0.1)
+  f_drift = NA, f_rd = 0.1,
+  scenario = FOCUS_Step_12_scenarios$names)
 {
+  if (n > 1 & is.na(i)) stop("Please specify the interval i if n > 1")
+
+  if (is.na(f_drift)) {
+    scenario = match.arg(scenario)
+    f_drift = FOCUS_Step_12_scenarios$drift[scenario, "1"] / 100
+    # For Step 2 we would/will select the reduced percentiles for multiple apps:
+    # if (n <= 8) {
+    #  f_drift = FOCUS_Step_12_scenarios$drift[scenario, as.character(n)] / 100
+    # } else {
+    #   f_drift = FOCUS_Step_12_scenarios$drift[scenario, ">8"] / 100
+    # }
+  }
+
   if (is.null(met)) {
+    cwsat = parent$cwsat
     mw_ratio = 1
     max_soil = 1
     max_ws = 1
     Koc = parent$Koc
     DT50_ws = parent$DT50_ws
   } else {
+    cwsat = met$cwsat
     mw_ratio = met$mw / parent$mw
     max_soil = met$max_soil
     max_ws = met$max_ws
@@ -41,17 +79,35 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
 
   # Rates for a single application
   eq_rate_drift_s = mw_ratio * max_ws * rate
+  # Parent only, or metabolite formed in soil:
   eq_rate_rd_s = mw_ratio * max_soil * rate
+  # Metabolite formed in water (this part is not documented in the Help files
+  # of the Steps 1/2 calculator):
+  if (!is.null(met)) {
+    eq_rate_rd_parent_s = mw_ratio * max_ws * rate
+    eq_rate_rd_s_tot = eq_rate_rd_s + eq_rate_rd_parent_s
+  } else {
+    eq_rate_rd_parent_s = NA
+    eq_rate_rd_s_tot = eq_rate_rd_s
+  }
 
   # Drift input
   input_drift_s = f_drift * eq_rate_drift_s / 10 # mg/m2
-  input_drift = n * input_drift_s
 
   # Runoff/drainage input
   ratio_field_wb = 10 # 10 m2 of field for each m2 of the waterbody
-  input_rd_s = f_rd * eq_rate_rd_s * ratio_field_wb / 10
+  input_rd_s = f_rd * eq_rate_rd_s_tot * ratio_field_wb / 10
   input_rd = n * input_rd_s
 
+  # No accumulation between multiple applications if 3 * DT50 < i
+  if (n > 1 && (3 * DT50_ws < i)) {
+    input_drift = input_drift_s
+    input_rd = input_rd_s
+  } else {
+    input_drift = n * input_drift_s
+    input_rd = n * input_rd_s
+  }
+
   # Fraction of compound entering the water phase via runoff/drainage
   depth_sw  = 0.3  # m
   depth_sed = 0.05 # m
@@ -85,25 +141,41 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
   PEC["1", "TWAECsw"] = (PEC_sw_0 + PEC["1", "PECsw"]) / 2
   PEC["1", "TWAECsed"] = (PEC_sed_0 + PEC["1", "PECsed"]) / 2
 
-  list(eq_rate_drift_s = eq_rate_drift_s,
-       eq_rate_rd_s = eq_rate_rd_s,
-       input_drift_s = input_drift_s,
-       input_rd_s = input_rd_s,
-       f_rd_sw = f_rd_sw, f_rd_sed = f_rd_sed,
-       PEC = PEC)
+  # Check if PEC_sw_max is above water solubility
+  PEC_sw_max = max(PEC[, "PECsw"])
+  if (PEC_sw_max > 1000 * cwsat) {
+    warning("The maximum PEC surface water exceeds the water solubility")
+  }
+
+  PEC_sed_max = max(PEC[, "PECsed"])
+
+  list(f_drift = f_drift,
+    eq_rate_drift_s = eq_rate_drift_s,
+    eq_rate_rd_s = eq_rate_rd_s,
+    eq_rate_rd_parent_s = eq_rate_rd_parent_s,
+    input_drift_s = input_drift_s,
+    input_rd_s = input_rd_s,
+    f_rd_sw = f_rd_sw, f_rd_sed = f_rd_sed,
+    PEC = PEC,
+    PEC_sw_max = PEC_sw_max,
+    PEC_sed_max = PEC_sed_max
+  )
 }
 
-#' Create an chemical compound object for FOCUS Step 1 and 2 calculations
+#' Create a chemical compound object for FOCUS Step 1 calculations
 #'
 #' @export
 #' @param cwsat Water solubility in mg/L
 #' @param DT50_ws Half-life in water/sediment systems in days
 #' @param Koc Partition coefficient between organic carbon and water
 #'   in L/kg.
+#' @param mw Molar weight in g/mol
+#' @param max_soil Maximum observed fraction (dimensionless) in soil
+#' @param max_ws Maximum observed fraction (dimensionless) in water/sediment
+#'   systems
 #' @return A list with the substance specific properties
-chent_focus_sw <- function(Koc, DT50_ws, cwsat = 1000)
+chent_focus_sw <- function(Koc, DT50_ws = NA, cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1)
 {
-  list(Koc = Koc, DT50_ws = DT50_ws, cwsat = cwsat)
+  list(Koc = Koc, DT50_ws = DT50_ws, cwsat = cwsat,
+       mw = mw, max_soil = max_soil, max_ws = max_ws)
 }
-
-
diff --git a/R/PELMO_runs.R b/R/PELMO_runs.R
index 0f134e0..b152b4d 100644
--- a/R/PELMO_runs.R
+++ b/R/PELMO_runs.R
@@ -3,7 +3,7 @@
 #' Per default, the runs are not only set up but also executed with FOCUS
 #' PELMO, the results are processed and returned. Currently, only FOCUS PELMO
 #' as installed on Linux (or other Unix systems)
-#' using the \code{\link{install_PELMO}} from the \code{PELMO.installeR} package
+#' using the \code{install_PELMO} from the \code{PELMO.installeR} package
 #' maintained on github is supported. In such installations, FOCUS PELMO is
 #' installed into the package installation directory of \code{PELMO.installeR}
 #' and run using \code{wine}.
@@ -33,6 +33,8 @@
 #' PELMO test results \url{http://esdac.jrc.ec.europa.eu/public_path/projects_data/focus/gw/models/pelmo/test-results/test_results_FOCUS_PELMO_5_5_3.doc}
 #' @export
 #' @examples
+#' # At the moment I can not run the examples, as my wine installation is not working
+#' \dontrun{
 #' # Reproduce the official test results for annual application of Pesticide D
 #' # to winter cereals at the day before emergence
 #' runs_1 <- list(
@@ -54,6 +56,7 @@
 #' PECgw_2 <- PELMO_runs(runs_2, psm_dir = system.file("testdata", package = "pfm"),
 #'   cores = 3, overwrite = TRUE)
 #' print(PECgw_2)
+#' }
 PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "auto",
                        execute = TRUE, cores = getOption("mc.cores", 2L),
                        evaluate = TRUE, overwrite = FALSE)
@@ -410,6 +413,9 @@ evaluate_PELMO <- function(runs, version = "5.5.3", PELMO_base = "auto")
 }
 
 #' Get the application interval in years from a psm file
+#'
+#' @param psm_file The path to the .psm file
+#' @param location_code The location code
 get_interval <- function(psm_file, location_code) {
   # How many years do we calculate (26, 46 or 66)?
   psm <- readLines(psm_file, encoding = "latin1")
@@ -433,6 +439,9 @@ get_interval <- function(psm_file, location_code) {
 }
 
 #' Sum up values according to FOCUS periods
+#' 
+#' @param annual The annual flux as obtained by \code{get_flux}
+#' @param interval The interval in years
 sum_periods <- function(annual, interval) {
   n_years <- switch(as.character(interval),
     "1" = 26,
@@ -446,6 +455,8 @@ sum_periods <- function(annual, interval) {
 }
 
 #' Get the flux of a chemical out of the FOCUS layer from a CHEM*.PLM file
+#' 
+#' @param chem_file The full path to a CHEM*.PLM file
 get_flux <- function(chem_file) {
   chem <- readLines(chem_file)
   lowest_focus_comp_lines <- grep("^    .     21      ", chem, value = TRUE)