Reorganise repository using standard package layout

14 files changed, 1242 insertions, 0 deletions

diff --git a/R/FOCUS_GW_scenarios_2012.R b/R/FOCUS_GW_scenarios_2012.R
new file mode 100644
index 0000000..49cf7dd
--- /dev/null
+++ b/R/FOCUS_GW_scenarios_2012.R
@@ -0,0 +1,11 @@
+#' A very small subset of the FOCUS Groundwater scenario defitions
+#' 
+#' Currently, only a small subset of the soil definitions are provided.
+#' 
+#' @name FOCUS_GW_scenarios_2012
+#' @references FOCUS (2012) Generic guidance for Tier 1 FOCUS ground water assessments. Version 2.1. 
+#'  FOrum for the Co-ordination of pesticde fate models and their USe. 
+#'  http://focus.jrc.ec.europa.eu/gw/docs/Generic_guidance_FOCV2_1.pdf
+#' @examples
+#' FOCUS_GW_scenarios_2012
+NULL
diff --git a/R/GUS.R b/R/GUS.R
new file mode 100644
index 0000000..4a6532d
--- /dev/null
+++ b/R/GUS.R
@@ -0,0 +1,98 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Groundwater ubiquity score based on Gustafson (1989)
+#' 
+#' The groundwater ubiquity score GUS is calculated according to
+#' the following equation
+#' \deqn{GUS = \log_{10} DT50_{soil} (4 - \log_{10} K_{oc})}{GUS = log10 DT50soil * (4 - log10 Koc)}
+#' 
+#' @references Gustafson, David I. (1989) Groundwater ubiquity score: a simple
+#' method for assessing pesticide leachability. \emph{Environmental
+#' toxicology and chemistry} \bold{8}(4) 339–57.
+#' @inheritParams endpoint
+#' @param DT50 Half-life of the chemical in soil. Should be a field
+#'   half-life according to Gustafson (1989). However, leaching to the sub-soil
+#'   can not completely be excluded in field dissipation experiments and Gustafson
+#'   did not refer to any normalisation procedure, but says the field study should
+#'   be conducted under use conditions.
+#' @param Koc The sorption constant normalised to organic carbon. Gustafson
+#'   does not mention the nonlinearity of the sorption constant commonly
+#'   found and usually described by Freundlich sorption, therefore it is 
+#'   unclear at which reference concentration the Koc should be observed
+#'   (and if the reference concentration would be in soil or in porewater).
+#' @param chent If a chent is given with appropriate information present in its
+#'   chyaml field, this information is used, with defaults specified below.
+#' @param degradation_value Which of the available degradation values should 
+#'   be used?
+#' @param lab_field Should laboratory or field half-lives be used? This
+#'   defaults to lab in this implementation, in order to avoid
+#'   double-accounting for mobility. If comparability with the original GUS
+#'   values given by Gustafson (1989) is desired, non-normalised first-order
+#'   field half-lives obtained under actual use conditions should be used.
+#' @param redox Aerobic or anaerobic degradation data
+#' @param sorption_value Which of the available sorption values should be used?
+#'   Defaults to Kfoc as this is what is generally available from the European
+#'   pesticide peer review process. These values generally use a reference
+#'   concentration of 1 mg/L in porewater, that means they would be expected to
+#'   be Koc values at a concentration of 1 mg/L in the water phase.
+#' @param degradation_aggregator Function for aggregating half-lives
+#' @param sorption_aggregator Function for aggregation Koc values
+#' @param ... Included in the generic to allow for further arguments later. Therefore
+#'   this also had to be added to the specific methods.
+#' @return A list with the DT50 and Koc used as well as the resulting score
+#'   of class GUS_result
+#' @author Johannes Ranke
+#' @export
+GUS <- function(...) UseMethod("GUS")
+
+#' @rdname GUS
+#' @export
+GUS.numeric <- function(DT50, Koc, ...) {
+  score <- log10(DT50) * (4 - log10(Koc))
+  res <- list(DT50 = DT50, Koc = Koc, score = score)
+  class(res) <- "GUS_result"
+  return(res)
+}
+
+#' @rdname GUS
+#' @export
+GUS.chent <- function(chent, 
+                      degradation_value = "DT50ref",
+                      lab_field = "laboratory",
+                      redox = "aerobic",
+                      sorption_value = "Kfoc", 
+                      degradation_aggregator = geomean,
+                      sorption_aggregator = geomean,
+                      ...)
+{
+  DT50 = soil_DT50(chent, lab_field = lab_field, redox = redox, 
+                   value = degradation_value, 
+                   aggregator = degradation_aggregator, signif = 5)
+  Koc = soil_Kfoc(chent, value = sorption_value, 
+                   aggregator = sorption_aggregator, signif = 5)
+  GUS.numeric(DT50, Koc)
+}
+
+#' @rdname GUS
+#' @export
+#' @param x An object of class GUS_result to be printed
+#' @param digits The number of digits used in the print method
+print.GUS_result = function(x, ..., digits = 1) {
+  cat("GUS: ", round(x$score, digits = 1), "\n")
+  cat("calculated from DT50 ", x$DT50, " and Koc ", x$Koc, "\n")
+}
diff --git a/R/PEC_soil.R b/R/PEC_soil.R
new file mode 100644
index 0000000..0263e47
--- /dev/null
+++ b/R/PEC_soil.R
@@ -0,0 +1,215 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+# Register global variables
+if(getRversion() >= '2.15.1') utils::globalVariables(c("destination", "study_type", "TP_identifier",
+                                                       "soil_scenario_data_EFSA_2015"))
+
+#' Calculate predicted environmental concentrations in soil
+#'
+#' This is a basic calculation of a contaminant concentration in bulk soil
+#' based on complete, instantaneous mixing. If an interval is given, an 
+#' attempt is made at calculating a long term maximum concentration using
+#' the concepts layed out for example in the PPR panel opinion (EFSA 2012).
+#' 
+#' This assumes that the complete load to soil during the time specified by
+#' 'interval' (typically 365 days) is dosed at once. As in the PPR panel
+#' opinion cited below (PPR panel 2012), only temperature correction using the
+#' Arrhenius equation is performed.  
+#' 
+#' Total soil and porewater PEC values for the scenarios as defined in the EFSA
+#' guidance (2015, p. 13) can easily be calculated.
+#' 
+#' @note If temperature information is available in the selected scenarios, as
+#'   e.g. in the EFSA scenarios, the DT50 for groundwater modelling
+#'   (destination 'PECgw') is taken from the chent object, otherwise the DT50
+#'   with destination 'PECsoil'.
+#' @importFrom methods is
+#' @param rate Application rate in units specified below
+#' @param rate_units Defaults to g/ha
+#' @param interception The fraction of the application rate that does not reach the soil
+#' @param mixing_depth Mixing depth in cm
+#' @param interval Period of the deeper mixing, defaults to 365, which is a year if
+#'   rate units are in days
+#' @param n_periods Number of periods to be considered for long term PEC calculations
+#' @param PEC_units Requested units for the calculated PEC. Only mg/kg currently supported
+#' @param PEC_pw_units Only mg/L currently supported
+#' @param tillage_depth Periodic (see interval) deeper mixing in cm
+#' @param chent An optional chent object holding substance specific information. Can 
+#'   also be a name for the substance as a character string
+#' @param DT50 If specified, overrides soil DT50 endpoints from a chent object
+#'   If DT50 is not specified here and not available from the chent object, zero
+#'   degradation is assumed
+#' @param Koc If specified, overrides Koc endpoints from a chent object
+#' @param Kom Calculated from Koc by default, but can explicitly be specified
+#'   as Kom here
+#' @param t_avg Averaging times for time weighted average concentrations
+#' @param scenarios If this is 'default', the DT50 will be used without correction
+#'   and soil properties as specified in the REACH guidance (R.16, Table
+#'   R.16-9) are used for porewater PEC calculations.  If this is "EFSA_2015",
+#'   the DT50 is taken to be a modelling half-life at 20°C and pF2 (for when
+#'   'chents' is specified, the DegT50 with destination 'PECgw' will be used),
+#'   and corrected using an Arrhenius activation energy of 65.4 kJ/mol. Also
+#'   model and scenario adjustment factors from the EFSA guidance are used.
+#' @param porewater Should equilibrium porewater concentrations be estimated
+#'   based on Kom and the organic carbon fraction of the soil instead of total
+#'   soil concentrations?  Based on equation (7) given in the PPR panel opinion
+#'   (EFSA 2012, p. 24) and the scenarios specified in the EFSA guidance (2015,
+#'   p. 13).
+#' @return The predicted concentration in soil
+#' @references EFSA Panel on Plant Protection Products and their Residues (2012)
+#'   Scientific Opinion on the science behind the guidance for scenario
+#'   selection and scenario parameterisation for predicting environmental
+#'   concentrations of plant protection products in soil. \emph{EFSA Journal}
+#'   \bold{10}(2) 2562, doi:10.2903/j.efsa.2012.2562
+#' 
+#'   EFSA (European Food Safety Authority) (2015) EFSA guidance document for
+#'   predicting environmental concentrations of active substances of plant
+#'   protection products and transformation products of these active substances
+#'   in soil. \emph{EFSA Journal} \bold{13}(4) 4093
+#'   doi:10.2903/j.efsa.2015.4093
+#' @author Johannes Ranke
+#' @export
+#' @examples
+#' PEC_soil(100, interception = 0.25)
+#' 
+#' # This is example 1 starting at p. 79 of the EFSA guidance (2015)
+#' PEC_soil(1000, interval = 365, DT50 = 250, t_avg = c(0, 21),
+#'                scenarios = "EFSA_2015")
+#' PEC_soil(1000, interval = 365, DT50 = 250, t_av = c(0, 21),
+#'                Kom = 1000, scenarios = "EFSA_2015", porewater = TRUE)
+#' 
+#' # The following is from example 4 starting at p. 85 of the EFSA guidance (2015)
+#' # Metabolite M2
+#' # Calculate total and porewater soil concentrations for tier 1 scenarios
+#' # Relative molar mass is 100/300, formation fraction is 0.7 * 1
+#' results_pfm <- PEC_soil(100/300 * 0.7 * 1 * 1000, interval = 365, DT50 = 250, t_avg = c(0, 21),
+#'                         scenarios = "EFSA_2015")
+#' results_pfm_pw <- PEC_soil(100/300 * 0.7 * 1000, interval = 365, DT50 = 250, t_av = c(0, 21),
+#'                            Kom = 100, scenarios = "EFSA_2015", porewater = TRUE)
+
+PEC_soil <- function(rate, rate_units = "g/ha", interception = 0,
+                     mixing_depth = 5, 
+                     PEC_units = "mg/kg", PEC_pw_units = "mg/L",
+                     interval = NA, n_periods = Inf,
+                     tillage_depth = 20, 
+                     chent = NA, 
+                     DT50 = NA, 
+                     Koc = NA, Kom = Koc / 1.724,
+                     t_avg = 0,
+                     scenarios = c("default", "EFSA_2015"),
+                     porewater = FALSE)
+{
+  rate_to_soil = (1 - interception) * rate
+  rate_units = match.arg(rate_units)
+  PEC_units = match.arg(PEC_units)
+  scenarios = match.arg(scenarios)
+  sce <- switch(scenarios, 
+    default = data.frame(rho = 1.5, T_arr = NA, theta_fc = 0.2, f_oc = 0.02, 
+                         f_sce = 1, f_mod = 1, row.names = "default"),
+    EFSA_2015 = if (porewater) soil_scenario_data_EFSA_2015[4:6, ] 
+                else soil_scenario_data_EFSA_2015[1:3, ]
+  )
+  n_sce = nrow(sce)
+
+  soil_volume = 100 * 100 * (mixing_depth/100)   # in m3
+  soil_mass = soil_volume * sce$rho * 1000  # in kg
+
+  # The following is C_T,ini from EFSA 2012, p. 22, but potentially with interception > 0
+  PEC_soil_ini = rate_to_soil * 1000 / soil_mass     # in mg/kg
+
+  # Decide which DT50 to take, or set degradation to zero if no DT50 available
+  if (is.na(DT50) & is(chent, "chent")) {
+    if (all(is.na(sce$T_arr))) {   # No temperature correction
+      DT50 <- subset(chent$soil_degradation_endpoints, destination == "PECsoil")$DT50
+    } else {
+      DT50 <- subset(chent$soil_degradation_endpoints, destination == "PECgw")$DT50
+    }
+    if (length(DT50) > 1) stop("More than one PECsoil DT50 in chent object")
+    if (length(DT50) == 0) DT50 <- Inf
+  }
+  k = log(2)/DT50
+
+  # Temperature correction of degradation (accumulation)
+  if (all(is.na(sce$T_arr))) {   # No temperature correction
+    f_T = 1
+  } else {
+    # Temperature correction as in EFSA 2012 p. 23
+    f_T = ifelse(sce$T_arr == 0, 
+                 0, 
+                 exp(- (65.4 / 0.008314) * (1/(sce$T_arr + 273.15) - 1/293.15)))
+  }
+
+  # X is the fraction left after one period (EFSA guidance p. 23)
+  X = exp(- k * f_T * interval)
+  
+  # f_accu is the fraction left after n periods (X + X^2 + ...)
+  f_accu = 0 
+  if (!is.na(interval)) {
+    if (n_periods == Inf) {
+      f_accu = X/(1 - X)
+    } else {
+      for (i in 1:n_periods) {
+        f_accu = f_accu + X^i
+      }
+    }
+  }
+
+  f_tillage = mixing_depth / tillage_depth
+
+  PEC_background = f_accu * f_tillage * PEC_soil_ini
+
+  PEC_soil = (1 + f_accu * f_tillage) * PEC_soil_ini
+
+  # Get porewater PEC if requested
+  if (porewater) {
+
+    # If Kom is not specified, try to get K(f)oc
+    if (is.na(Kom)) {
+      # If Koc not specified, try to get K(f)oc from chent
+      if (is.na(Koc) & is(chent, "chent")) {
+        Koc <- soil_Kfoc(chent)
+      } 
+      Kom <- Koc / 1.724
+    }
+
+    if (is.na(Kom)) stop("No Kom information specified")
+
+    PEC_soil = PEC_soil/((sce$theta_fc/sce$rho) + sce$f_om * Kom)
+  }
+
+  # Scenario adjustment factors
+  PEC_soil_sce = PEC_soil * sce$f_sce
+
+  # Model adjustment factors
+  PEC_soil_sce_mod = PEC_soil_sce * sce$f_mod
+
+  result <- matrix(NA, ncol = n_sce, nrow = length(t_avg),
+                   dimnames = list(t_avg = t_avg, scenario = rownames(sce)))
+  
+  result[1, ] <- PEC_soil_sce_mod
+
+  for (i in seq_along(t_avg)) {
+    t_av_i <- t_avg[i]
+    if (t_av_i > 0) {
+      # Equation 10 from p. 24 (EFSA 2015)
+      result[i, ] <- PEC_soil_sce_mod/(t_av_i * f_T * k) * (1 - exp(- f_T * k * t_av_i))
+    }
+  }
+
+  return(result)
+}
diff --git a/R/PEC_sw_drainage_UK.R b/R/PEC_sw_drainage_UK.R
new file mode 100644
index 0000000..e53f179
--- /dev/null
+++ b/R/PEC_sw_drainage_UK.R
@@ -0,0 +1,66 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Calculate initial predicted environmental concentrations in surface water due to drainage using the UK method
+#'
+#' This implements the method specified in the UK data requirements handbook and was checked against the spreadsheet
+#' published on the CRC website
+#'
+#' @param rate Application rate in g/ha
+#' @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 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 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
+#' @return The predicted concentration in surface water in µg/L
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' PEC_sw_drainage_UK(150, Koc = 100)
+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
+
+  if (!missing(latest_application)) {
+    lct <- Sys.getlocale("LC_TIME")
+    tmp <- Sys.setlocale("LC_TIME", "C")
+    latest <- as.Date(paste(latest_application, "1999"), "%d %b %Y")
+    tmp <- Sys.setlocale("LC_TIME", lct)
+    degradation_time <- as.numeric(difftime(as.Date("1999-10-01"), units = "days", latest))
+    if (!missing(soil_DT50)) {
+      k = log(2)/soil_DT50
+      as.Date(paste(latest_application, "1999"), "%d %B %Y")
+
+      amount_available <- amount_available * exp(-k * degradation_time)
+      if (!missing(model)) stop("You already supplied a soil_DT50 value, implying SFO kinetics")
+    }
+    if (!missing(model)) {
+      fraction_left <- pfm_degradation(model, parms = model_parms, 
+                                       times = degradation_time)[1, "parent"]
+      amount_available <- fraction_left * amount_available
+    }
+  } 
+
+  volume = 130000 # L/ha
+  PEC = 1e6 * (percentage_lost/100) * amount_available / volume
+  return(PEC)
+}
diff --git a/R/PEC_sw_drift.R b/R/PEC_sw_drift.R
new file mode 100644
index 0000000..261706c
--- /dev/null
+++ b/R/PEC_sw_drift.R
@@ -0,0 +1,65 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Calculate predicted environmental concentrations in surface water due to drift
+#'
+#' This is a basic, vectorised form of a simple calculation of a contaminant
+#' concentration in surface water based on complete, instantaneous mixing
+#' with input via spray drift.
+#'
+#' @param rate Application rate in units specified below
+#' @param applications Number of applications for selection of drift percentile
+#' @param drift_percentages Percentage drift values for which to calculate PECsw.
+#'   'drift_data' and 'distances' if not NULL.
+#' @param drift_data Source of drift percentage data
+#' @param crop Crop name (use German names for JKI data), defaults to "Ackerbau"
+#' @param distances The distances in m for which to get PEC values
+#' @param water_depth Depth of the water body in cm
+#' @param rate_units Defaults to g/ha
+#' @param PEC_units Requested units for the calculated PEC. Only µg/L currently supported
+#' @return The predicted concentration in surface water
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' PEC_sw_drift(100)
+PEC_sw_drift <- function(rate, 
+                         applications = 1,
+                         water_depth = 30, 
+                         drift_percentages = NULL,
+                         drift_data = "JKI",
+                         crop = "Ackerbau",
+                         distances = c(1, 5, 10, 20),
+                         rate_units = "g/ha",
+                         PEC_units = "\u00B5g/L")
+{
+  rate_units <- match.arg(rate_units)
+  PEC_units <- match.arg(PEC_units)
+  drift_data <- match.arg(drift_data)
+  water_volume <- 100 * 100 * (water_depth/100) * 1000   # in L (for 1 ha)
+  PEC_sw_overspray <- rate * 1e6 / water_volume          # in µg/L
+  dist_index <- as.character(distances)
+
+  if (is.null(drift_percentages)) {
+    drift_percentages <- pfm::drift_data_JKI[[applications]][dist_index, crop]
+    names(drift_percentages) <- paste(dist_index, "m")
+  } else {
+    names(drift_percentages) <- paste(drift_percentages, "%")
+  }
+
+  PEC_sw_drift <- PEC_sw_overspray * drift_percentages / 100
+  return(PEC_sw_drift)
+}
diff --git a/R/PEC_sw_sed.R b/R/PEC_sw_sed.R
new file mode 100644
index 0000000..a4a3d6a
--- /dev/null
+++ b/R/PEC_sw_sed.R
@@ -0,0 +1,50 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Calculate predicted environmental concentrations in sediment from surface
+#' water concentrations
+#'
+#' The method 'percentage' is equivalent to what is used in the CRD spreadsheet
+#' PEC calculator
+#'
+#' @param PEC_sw Numeric vector or matrix of surface water concentrations in µg/L for
+#'   which the corresponding sediment concentration is to be estimated
+#' @param percentage The percentage in sediment, used for the percentage method
+#' @param method The method used for the calculation
+#' @param sediment_depth Depth of the sediment layer
+#' @param water_depth Depth of the water body in cm
+#' @param sediment_density The density of the sediment in L/kg (equivalent to
+#'   g/cm3)
+#' @param PEC_sed_units The units of the estimated sediment PEC value
+#' @return The predicted concentration in sediment
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' PEC_sw_sed(PEC_sw_drift(100, distances = 1), percentage = 50)
+PEC_sw_sed <- function(PEC_sw, percentage = 100, method = "percentage", 
+                       sediment_depth = 5, water_depth = 30,
+                       sediment_density = 1.3,
+                       PEC_sed_units = c("\u00B5g/kg", "mg/kg"))
+{
+  method = match.arg(method)
+  PEC_sed_units = match.arg(PEC_sed_units)
+  if (method == "percentage") {
+    PEC_sed = PEC_sw * (percentage/100) * (water_depth / sediment_depth) * (1 / sediment_density)
+    if (PEC_sed_units == "mg/kg") PEC_sed <- PEC_sed / 1000
+  }
+  return(PEC_sed)
+}
diff --git a/R/SFO_actual_twa.R b/R/SFO_actual_twa.R
new file mode 100644
index 0000000..7facb6a
--- /dev/null
+++ b/R/SFO_actual_twa.R
@@ -0,0 +1,36 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Actual and maximum moving window time average concentrations for SFO kinetics
+#'
+#' @param DT50 The half-life.
+#' @param times The output times, and window sizes for time weighted average concentrations
+#' @export
+#' @author Johannes Ranke
+#' @source FOCUS (2014) Generic Guidance for Estimating Persistence and Degradation
+#'   Kinetics from Environmental Fate Studies on Pesticides in EU Registratin, Version 1.1, 
+#'   18 December 2014, p. 251
+#' @examples
+#' SFO_actual_twa(10)
+SFO_actual_twa <- function(DT50 = 1000, times = c(0, 1, 2, 4, 7, 14, 21, 28, 42, 50, 100))
+{
+  k = log(2)/DT50
+  result <- data.frame(actual = 1 * exp(-k * times),
+                       twa = (1 - exp(-k * times))/(k * times),
+                       row.names = times)
+  return(result)
+}
diff --git a/R/SSLRC_mobility_classification.R b/R/SSLRC_mobility_classification.R
new file mode 100644
index 0000000..deda5cf
--- /dev/null
+++ b/R/SSLRC_mobility_classification.R
@@ -0,0 +1,42 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Determine the SSLRC mobility classification for a chemical substance from its Koc
+#'
+#' This implements the method specified in the UK data requirements handbook and was 
+#' checked against the spreadsheet published on the CRC website
+#'
+#' @param Koc The sorption coefficient normalised to organic carbon in L/kg
+#' @return A list containing the classification and the percentage of the
+#'   compound transported per 10 mm drain water
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' SSLRC_mobility_classification(100)
+SSLRC_mobility_classification <- function(Koc)
+{
+  if (!is.numeric(Koc) | length(Koc) != 1) stop("Please give a single number")
+  result <- list("Non mobile", 0.01)
+  if (Koc < 4000) result <- list("Slightly mobile", 0.02)
+  if (Koc < 1000) result <- list("Slightly mobile", 0.5)
+  if (Koc < 500) result <- list("Moderately mobile", 0.7)
+  if (Koc < 75) result <- list("Mobile", 1.9)
+  if (Koc < 15) result <- list("Very mobile", 1.9)
+  names(result) <- c("Mobility classification", 
+                     "Percentage drained per mm of drain water")
+  return(result)
+}
diff --git a/R/TOXSWA_cwa.R b/R/TOXSWA_cwa.R
new file mode 100644
index 0000000..c5ddce9
--- /dev/null
+++ b/R/TOXSWA_cwa.R
@@ -0,0 +1,380 @@
+# Copyright (C) 2014,2015,2016  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Read TOXSWA surface water concentrations
+#'
+#' Read TOXSWA hourly concentrations of a chemical substance in a specific
+#' segment of a TOXSWA surface water body. Per default, the data for the last
+#' segment are imported. As TOXSWA 4 reports the values at the end of the hour
+#' (ConLiqWatLayCur) in its summary file, we use this value as well instead
+#' of the hourly averages (ConLiqWatLay).
+#'
+#' @param filename The filename of the cwa file (TOXSWA 2.x.y  or similar) or the
+#'  out file (FOCUS TOXSWA 4, i.e. TOXSWA 4.4.2 or similar).
+#' @param basedir The path to the directory where the cwa file resides.
+#' @param zipfile Optional path to a zip file containing the cwa file.
+#' @param segment The segment for which the data should be read. Either "last", or 
+#'   the segment number.
+#' @param total Set this to TRUE in order to read total concentrations as well. This is 
+#'   only necessary for .out files as generated by TOXSWA 4.4.2 or similar, not for .cwa
+#'   files. For .cwa files, the total concentration is always read as well.
+#' @param substance For TOXSWA 4 .out files, the default value "parent" leads
+#'   to reading concentrations of the parent compound. Alternatively, the substance
+#'   of interested can be selected by its code name.
+#' @param windows Numeric vector of width of moving windows in days, for calculating
+#'   maximum time weighted average concentrations and areas under the curve.
+#' @param thresholds Numeric vector of threshold concentrations in µg/L for
+#'   generating event statistics.  
+#' @importFrom readr read_fwf fwf_empty
+#' @return An instance of an R6 object of class
+#' \code{\link{TOXSWA_cwa}}.
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' H_sw_D4_pond  <- read.TOXSWA_cwa("00001p_pa.cwa",
+#'                                  basedir = "SwashProjects/project_H_sw/TOXSWA",
+#'                                  zipfile = system.file("testdata/SwashProjects.zip",
+#'                                                        package = "pfm"))
+read.TOXSWA_cwa <- function(filename, basedir = ".", zipfile = NULL, 
+                            segment = "last", substance = "parent",
+                            total = FALSE,
+                            windows = NULL, thresholds = NULL)
+{
+  if (!missing(filename)) {
+    cwa <- TOXSWA_cwa$new(filename, basedir, zipfile, 
+                          substance = substance, total = total)
+    if (!is.null(windows[1])) cwa$moving_windows(windows)
+    if (!is.null(thresholds[1])) cwa$get_events(thresholds)
+    invisible(cwa)
+  } else {
+    message("You need to specify a filename for the cwa file to be read")
+  }
+}
+
+#' Plot TOXSWA surface water concentrations
+#'
+#' Plot TOXSWA hourly concentrations of a chemical substance in a specific
+#' segment of a TOXSWA surface water body.
+#'
+#' @import graphics
+#' @param x The TOXSWA_cwa object to be plotted.
+#' @param xlab,ylab Labels for x and y axis.
+#' @param time_column What should be used for the time axis. If "t_firstjan" is chosen,
+#'   the time is given in days relative to the first of January in the first year.
+#' @param add Should we add to an existing plot?
+#' @param total Should the total concentration in water be plotted, including substance sorbed
+#'   to suspended matter?
+#' @param LC_TIME Specification of the locale used to format dates
+#' @param ... Further arguments passed to \code{plot} if we are not adding to an existing plot
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' H_sw_D4_pond  <- read.TOXSWA_cwa("00001p_pa.cwa",
+#'                                  basedir = "SwashProjects/project_H_sw/TOXSWA",
+#'                                  zipfile = system.file("testdata/SwashProjects.zip",
+#'                                              package = "pfm"))
+#' plot(H_sw_D4_pond)
+plot.TOXSWA_cwa <- function(x, time_column = c("datetime", "t", "t_firstjan", "t_rel_to_max"),
+                            xlab = "default", ylab = "default",
+                            add = FALSE,
+                            total = FALSE, LC_TIME = "C", ...)
+{
+  time_column = match.arg(time_column)
+  cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L")
+  lct <- Sys.getlocale("LC_TIME")
+  tmp <- Sys.setlocale("LC_TIME", LC_TIME)
+  if (xlab == "default") {
+    xlab = switch(time_column,
+                datetime = "Time",
+                t = "Time [days]",
+                t_firstjan = "Time since first of January [days]",
+                t_rel_to_max = "Time relative to maximum concentration [days]")  
+  }
+  if (ylab == "default") {
+    ylab = paste( ifelse(total, "Total concentration", "Concentration"), "[\u03bcg/L]")
+  }
+  if (add) {
+    lines(x$cwas[c(time_column, cwa_column)], xlab = xlab, ylab = ylab, ...)
+  } else{
+    if (time_column == "datetime") {
+      plot(x$cwas$datetime, x$cwas[[cwa_column]], type = "l",
+           xlab = xlab, ylab = ylab, xaxt = "n", ...)
+      seq_x <- seq(min(x$cwas$datetime), max(x$cwas$datetime), by = "quarter")
+      axis(1, at = seq_x, labels = format(seq_x, "%b %Y"))
+    } else {
+      plot(x$cwas[c(time_column, cwa_column)], type = "l",
+           xlab = xlab, ylab = ylab, ...)
+    }
+  }
+  tmp <- Sys.setlocale("LC_TIME", lct)
+}
+
+#' R6 class for holding TOXSWA cwa concentration data and associated statistics
+#'
+#' An R6 class for holding TOXSWA cwa concentration data and some associated statistics.
+#' Usually, an instance of this class will be generated by \code{\link{read.TOXSWA_cwa}}.
+#'
+#' @docType class
+#' @importFrom R6 R6Class
+#' @export
+#' @format An \code{\link{R6Class}} generator object.
+#' @field filename Length one character vector.
+#' @field basedir Length one character vector.
+#' @field segment Length one integer, specifying for which segment the cwa data were read.
+#' @field cwas Dataframe holding the concentrations.
+#' @field events List of dataframes holding the event statistics for each threshold.
+#' @field windows Matrix of maximum time weighted average concentrations (TWAC_max)
+#'   and areas under the curve in µg/day * h (AUC_max_h) or µg/day * d (AUC_max_d) 
+#'   for the requested moving window sizes in days.
+#' @section Methods:
+#' \describe{
+#'   \item{\code{get_events(threshold, total = FALSE)}}{
+#'         Populate a datataframe with event information for the specified threshold value
+#'         in µg/L. If \code{total = TRUE}, the total concentration including the amount
+#'         adsorbed to suspended matter will be used. The resulting dataframe is stored in the 
+#'         \code{events} field of the object.
+#'         }
+#'   \item{\code{moving_windows(windows, total = FALSE)}}{
+#'         Add to the \code{windows} field described above.
+#'         Again, if \code{total = TRUE}, the total concentration including the amount
+#'         adsorbed to suspended matter will be used.
+#'         }
+#' }
+#' @examples
+#' H_sw_R1_stream  <- read.TOXSWA_cwa("00003s_pa.cwa",
+#'                                  basedir = "SwashProjects/project_H_sw/TOXSWA",
+#'                                  zipfile = system.file("testdata/SwashProjects.zip",
+#'                                              package = "pfm"))
+#' H_sw_R1_stream$get_events(c(2, 10))
+#' H_sw_R1_stream$moving_windows(c(7, 21))
+#' print(H_sw_R1_stream)
+#' @keywords data
+
+TOXSWA_cwa <- R6Class("TOXSWA_cwa",
+  public = list(
+    filename = NULL,
+    basedir = NULL,
+    zipfile = NULL,
+    segment = NULL,
+    substance = NULL,
+    cwas = NULL,
+    windows = NULL,
+    events = list(),
+    initialize = function(filename, basedir, zipfile = NULL, 
+                          segment = "last", substance = "parent", total = FALSE) {
+      self$filename <- filename
+      self$basedir <- basedir
+      self$zipfile <- zipfile
+      if (!is.null(zipfile)) {
+        try(file_connection <- unz(zipfile, paste0(basedir, "/", filename)))
+      } else {
+        try(file_connection <- file(file.path(basedir, filename), "rt"))
+      }
+      
+
+      if (grepl(".cwa$", filename)) {
+        # cwa file from FOCUS TOXSWA 3 (TOXSWA 2.x.y)
+        cwa_all_segments <- try(
+          read.table(file_connection,
+          sep = "", skip = 40,
+          encoding = "UTF-8",
+          colClasses = c("character", "numeric",
+                         "integer", rep("numeric", 5)),
+          col.names = c("datetime", "t", "segment",
+                        "xcd", "cwa_tot", "cwa", "Xss", "Xmp")))
+
+        if (is.null(zipfile)) close(file_connection) # only needed for files
+
+        if (!inherits(cwa_all_segments, "try-error")) {
+          available_segments = 1:max(cwa_all_segments$segment)
+          if (segment == "last") segment = max(available_segments)
+          if (!segment %in% available_segments) stop("Invalid segment specified")
+          self$segment <- segment
+          cwa <- subset(cwa_all_segments, segment == self$segment, 
+                        c("datetime", "t", "segment", "cwa", "cwa_tot"))
+          lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C")
+          cwa$datetime <- strptime(cwa$datetime, "%d-%b-%Y-%H:%M")
+          Sys.setlocale("LC_TIME", lct)
+          startyear = format(cwa$datetime[1], "%Y") 
+          firstjan <- strptime(paste0(startyear, "-01-01"), "%Y-%m-%d")
+          cwa$t_firstjan <- as.numeric(difftime(cwa$datetime, 
+                                                      firstjan, units = "days"))
+
+          t_max = cwa[which.max(cwa$cwa), "t"]
+          cwa$t_rel_to_max = cwa$t - t_max
+          cwa$cwa_mug_per_L <- cwa$cwa * 1000
+          cwa$cwa_tot_mug_per_L <- cwa$cwa_tot * 1000
+          self$cwas <- cwa[c("datetime", "t", "t_firstjan",
+                                          "t_rel_to_max",
+                                          "cwa_mug_per_L", 
+                                          "cwa_tot_mug_per_L")]
+        } else {
+          stop("Could not read ", filename)
+        }
+      } else {
+        # out file from FOCUS TOXSWA 4 (TOXSWA 4.4.2 or similar)
+        outfile <- try(readLines(file_connection))
+
+        close(file_connection) # only needed for files
+
+        if (inherits(outfile, "try-error")) {
+          stop("Could not read ", filename)
+        } else {
+          # Get the substance name(s)
+          sub_lines <- grep(".*0.000.*ConLiqWatLayCur_", outfile[1:50], value = TRUE)
+          substances <- gsub(".*ConLiqWatLayCur_(.*?) *[0-9].*", "\\1", sub_lines)
+          if (!substance %in% c("parent", substances)) {
+            stop("No data for substance ", substance, " present in the .out file.")
+          }
+          
+          # Generate field name for the concentrations at end of hour for the
+          # substance of interest
+          if (substance == "parent") {
+            cwa_string = paste0("ConLiqWatLayCur_", substances[1])
+          } else {
+            cwa_string = paste0("ConLiqWatLayCur_", substance)
+          }
+
+          cwa_lines <- grep(cwa_string, outfile, value = TRUE) 
+
+          cwa_all_segments <- read_fwf(paste(cwa_lines, collapse = "\n"), 
+                                       fwf_empty(paste(tail(cwa_lines), collapse = "\n")))
+
+          # Append time "-00h00" to datetime in first row, as this is not (always?) present
+          # in the line ConLiqWatLayCur
+          if (nchar(cwa_all_segments[1, "X2"]) == 11) {
+            cwa_all_segments[1, "X2"] = paste0(cwa_all_segments[1, "X2"], "-00h00")
+          }
+
+          available_segments = 1:(ncol(cwa_all_segments) - 3)
+          if (segment == "last") segment = max(available_segments)
+          if (!segment %in% available_segments) stop("Invalid segment specified")
+          self$segment <- segment
+          cwa <- data.frame(
+            datetime = as.character(cwa_all_segments$X2),
+            t = cwa_all_segments$X1,
+            cwa = cwa_all_segments[[3 + segment]]
+          )
+          if (total) {
+            cwa_tot_lines <- outfile[grep("ConSysWatLay_", outfile)] # hourly total conc.
+            cwa_tot_all_segments <- read.table(text = cwa_lines) 
+            cwa$cwa_tot = cwa_tot_all_segments[[3 + segment]]
+          }
+          lct <- Sys.getlocale("LC_TIME"); Sys.setlocale("LC_TIME", "C")
+          cwa$datetime <- strptime(cwa$datetime, "%d-%b-%Y-%Hh%M")
+          Sys.setlocale("LC_TIME", lct)
+
+          startyear = format(cwa$datetime[1], "%Y") 
+          firstjan <- strptime(paste0(startyear, "-01-01"), "%Y-%m-%d")
+          cwa$t_firstjan <- as.numeric(difftime(cwa$datetime, 
+                                                firstjan, units = "days"))
+          t_max = cwa[which.max(cwa$cwa), "t"]
+          cwa$t_rel_to_max = cwa$t - t_max
+          cwa$cwa_mug_per_L <- cwa$cwa * 1000
+          self$cwas <- cwa[c("datetime", "t", "t_firstjan",
+                                          "t_rel_to_max",
+                                          "cwa_mug_per_L")]
+          if (total) {
+            self$cwas$cwa_tot_mug_per_L <- cwa$cwa_tot * 1000
+          }
+        }
+      }
+    },
+    moving_windows = function(windows, total = FALSE) {
+      window_names = paste(windows, "days")
+      n = length(window_names)
+      self$windows <- data.frame(window = window_names,
+                                 max_TWAC = numeric(n),
+                                 max_AUC_h = numeric(n),
+                                 max_AUC_d = numeric(n))
+      if (missing(windows)) {
+        stop("You need to specify at least one moving window size in days")
+      }
+      cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L")
+      for (i in seq_along(windows)) {
+        window_size = windows[i]
+        filter_size = window_size * 24
+        max_TWAC = max(filter(self$cwas[cwa_column], rep(1/filter_size, filter_size), 
+                            "convolution"), na.rm = TRUE)
+        max_AUC_h = max_TWAC * filter_size
+        max_AUC_d = max_TWAC * window_size
+        
+        self$windows[i, -1] = c(max_TWAC, max_AUC_h, max_AUC_d)
+      }
+      invisible(self)
+    },
+    get_events = function(thresholds, total = FALSE) {
+      if (missing(thresholds)) {
+        stop("You need to specify at least one threshold concentration in \u03bcg/L")
+      }
+      for (threshold in thresholds) {
+        events = data.frame(t_start = numeric(), cwa_max = numeric(), 
+                            duration = numeric(), pre_interval = numeric(), 
+                            AUC_h = numeric(), AUC_d = numeric())
+        cwa_column = ifelse(total, "cwa_tot_mug_per_L", "cwa_mug_per_L")
+
+        event_end = 0
+        event = FALSE
+        event_max = 0
+        event_nr = 0
+        n_rows = nrow(self$cwas)
+        for (i in 1:n_rows) {
+          cwa_cur = self$cwas[i, cwa_column]
+          if (event == FALSE) {
+            if (cwa_cur > threshold) {
+              event_start = self$cwas[i, "t"]
+              pre_interval = event_start - event_end
+              i_start = i
+              event = TRUE
+              event_max = cwa_cur
+            }
+          } else {
+            if (cwa_cur > event_max) event_max = cwa_cur
+            if (cwa_cur < threshold || i == n_rows) {
+              event_nr = event_nr + 1
+              i_end = i
+              if (i == n_rows) i_end = i
+              event_end = self$cwas[i_end, "t"]
+              event_length = event_end - event_start
+              event_cwas <- self$cwas[i_start:(i_end - 1), cwa_column]
+              event_AUC_h = sum(event_cwas)
+              event_AUC_d = event_AUC_h / 24
+              events[event_nr, ] = c(event_start, event_max, event_length,
+                                     pre_interval, event_AUC_h, event_AUC_d)
+              event = FALSE
+            } 
+          }
+        }
+        
+        self$events[[as.character(threshold)]] <- events
+      }
+      invisible(self)
+    },
+    print = function() {
+      cat("<TOXSWA_cwa> data from file", self$filename, "segment", self$segment, "\n")
+      print(head(self$cwas))
+      cat("Moving window analysis\n")
+      print(self$windows)
+      for (threshold in names(self$events)) {
+        cat("Event statistics for threshold", threshold, "\n")
+        if (nrow(self$events[[threshold]]) == 0) cat("No events found\n")
+        else print(self$events[[threshold]])
+      }
+    }
+  )
+)
+# vim: set ts=2 sw=2 expandtab:
diff --git a/R/drift_data_JKI.R b/R/drift_data_JKI.R
new file mode 100644
index 0000000..888d3b5
--- /dev/null
+++ b/R/drift_data_JKI.R
@@ -0,0 +1,44 @@
+#' Deposition from spray drift expressed as percent of the applied dose as
+#' published by the JKI
+#' 
+#' Deposition from spray drift expressed as percent of the applied dose as
+#' published by the German Julius-Kühn Institute (JKI).
+#' 
+#' The data were extracted from the spreadsheet cited below using the R code
+#' given in the example section. The spreadsheet is not included in the package
+#' as its licence is not clear.
+#' 
+#' 
+#' @name drift_data_JKI
+#' @docType data
+#' @format A list currently containing matrices with spray drift percentage
+#' data for field crops (Ackerbau), and Pome/stone fruit, early and late
+#' (Obstbau frueh, spaet).
+#' @source JKI (2010) Spreadsheet 'Tabelle der Abdrifteckwerte.xls', retrieved
+#' from
+#' http://www.jki.bund.de/no_cache/de/startseite/institute/anwendungstechnik/abdrift-eckwerte.html
+#' on 2015-06-11
+#' @keywords datasets
+#' @examples
+#' 
+#' \dontrun{
+#'   # This is the code that was used to extract the data
+#'   library(readxl)
+#'   abdrift_path <- "inst/extdata/Tabelle der Abdrifteckwerte.xls"
+#'   JKI_crops <- c("Ackerbau", "Obstbau frueh", "Obstbau spaet")
+#'   names(JKI_crops) <- c("Field crops", "Pome/stone fruit, early", "Pome/stone fruit, late")
+#'   drift_data_JKI <- list()
+#' 
+#'   for (n in 1:8) {
+#'     drift_data_raw <- read_excel(abdrift_path, sheet = n + 1, skip = 2)
+#'     drift_data <- as.matrix(drift_data_raw[1:9, 2:4]) 
+#'     dimnames(drift_data) <- list(distance = as.integer(drift_data_raw[1:9, 1]),
+#'                                             crop = JKI_crops)
+#'     drift_data_JKI[[n]] <- drift_data
+#'   }
+#'   save(drift_data_JKI, file = "data/drift_data_JKI.RData")
+#' }
+#' 
+#' # And this is the resulting data
+#' drift_data_JKI
+NULL
diff --git a/R/endpoint.R b/R/endpoint.R
new file mode 100644
index 0000000..56a314b
--- /dev/null
+++ b/R/endpoint.R
@@ -0,0 +1,109 @@
+#' Retrieve endpoint information from the chyaml field of a chent object
+#'
+#' R6 class objects of class chent represent chemical entities
+#' and can hold a list of information loaded from a chemical yaml file in their
+#' chyaml field. Such information is extracted and optionally aggregated by
+#' this function.
+#' 
+#' The functions \code{soil_*} are functions to extract soil specific endpoints.
+#' For the Freundlich exponent, the capital letter \code{N} is used in order to
+#' facilitate dealing with such data in R. In pesticide fate modelling, this
+#' exponent is often called 1/n.
+#'
+#' @export
+#' @param chent The chent object to get the information from
+#' @param medium The medium for which information is sought
+#' @param type The information type
+#' @param lab_field If not NA, do we want laboratory or field endpoints
+#' @param redox If not NA, are we looking for aerobic or anaerobic data
+#' @param value The name of the value we want. The list given in the 
+#'   usage section is not exclusive
+#' @param aggregator The aggregator function. Can be mean, 
+#'   \code{\link{geomean}}, or identity, for example.
+#' @param raw Should the number(s) be returned as stored in the chent
+#'   object (could be a character value) to retain original information
+#'   about precision?
+#' @param signif How many significant digits do we want
+#' @return The result from applying the aggregator function to
+#'   the values converted to a numeric vector, rounded to the
+#'   given number of significant digits, or, if raw = TRUE,
+#'   the values as a character value, retaining any implicit
+#'   information on precision that may be present.
+#' 
+endpoint <- function(chent, 
+                     medium = "soil",
+                     type = c("degradation", "sorption"), 
+                     lab_field = c(NA, "laboratory", "field"),
+                     redox = c(NA, "aerobic", "anaerobic"),
+                     value = c("DT50ref", "Kfoc", "N"),
+                     aggregator = geomean,
+                     raw = FALSE,
+                     signif = 3)
+{
+  if (!is(chent, "chent")) {
+    stop("Please supply a chent object as created using the package 'chents' available from jrwb.de")
+  }
+  ep_list <- chent$chyaml[[medium]][[type]] 
+  if (!is.na(lab_field[1])) {
+    ep_list <- ep_list[[lab_field]]
+  }
+  if (!is.na(redox[1])) {
+    ep_list <- ep_list[[redox]]
+  }
+  values <- ep_list$data[[value]]
+  if (raw) return(values)
+  else return(signif(aggregator(as.numeric(values)), signif))
+}
+
+#' @inheritParams endpoint
+#' @rdname endpoint
+#' @export
+soil_DT50 <- function(chent, aggregator = geomean, signif = 3, 
+                      lab_field = "laboratory", value = "DT50ref",
+                      redox = "aerobic", raw = FALSE) {
+  ep <- endpoint(chent, medium = "soil", type = "degradation", 
+                 lab_field = "laboratory", redox = redox,
+                 signif = signif,
+                 value = value, aggregator = aggregator, raw = raw)
+  return(ep)
+}
+
+#' @inheritParams endpoint
+#' @rdname endpoint
+#' @export
+soil_Kfoc <- function(chent, aggregator = geomean, signif = 3, 
+                      value = "Kfoc", raw = FALSE) {
+  ep <- endpoint(chent, medium = "soil", type = "sorption", 
+                 signif = signif,
+                 value = value, aggregator = aggregator, raw = raw)
+  return(ep)
+}
+
+#' @inheritParams endpoint
+#' @rdname endpoint
+#' @export
+soil_N <- function(chent, aggregator = mean, signif = 3, raw = FALSE) {
+  ep <- endpoint(chent, medium = "soil", type = "sorption", 
+                 signif = signif,
+                 value = "N", aggregator = aggregator, raw = raw)
+  return(ep)
+}
+
+#' @inheritParams endpoint
+#' @rdname endpoint
+#' @param values The values to be returned
+#' @param aggregators A named vector of aggregator functions to be used
+#' @export
+soil_sorption <- function(chent, values = c("Kfoc", "N"), 
+                          aggregators = c(Kfoc = geomean, Koc = geomean, N = mean), 
+                          signif = c(Kfoc = 3, N = 3),
+                          raw = FALSE) {
+  res <- sapply(values, 
+                function(x) {
+                  endpoint(chent, medium = "soil", type = "sorption",
+                           signif = signif[[x]],
+                           value = x, aggregator = aggregators[[x]], raw = raw)
+                }
+  )
+  return(res)
+}
diff --git a/R/geomean.R b/R/geomean.R
new file mode 100644
index 0000000..626829b
--- /dev/null
+++ b/R/geomean.R
@@ -0,0 +1,38 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Calculate the geometric mean
+#'
+#' Based on some posts in a thread on Stackoverflow
+#' \url{http://stackoverflow.com/questions/2602583/geometric-mean-is-there-a-built-in}
+#' This function checks for negative values, removes NA values per default and
+#' returns 0 if at least one element of the vector is 0.
+#'
+#' @param x Vector of numbers
+#' @param na.rm Should NA values be omitted?
+#' @return The geometric mean
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' geomean(c(1, 3, 9))
+#' geomean(c(1, 3, NA, 9))
+#' \dontrun{geomean(c(1, -3, 9)) # returns an error}
+geomean = function(x, na.rm = TRUE){
+  if (any(is.na(x)) & na.rm == FALSE) stop("Removal of NA values was not requested")
+  if (any(x < 0, na.rm = na.rm)) stop("Only defined for positive numbers")
+  exp(mean(log(x), na.rm = na.rm))
+}
diff --git a/R/pfm_degradation.R b/R/pfm_degradation.R
new file mode 100644
index 0000000..6c8610e
--- /dev/null
+++ b/R/pfm_degradation.R
@@ -0,0 +1,48 @@
+# Copyright (C) 2015  Johannes Ranke
+# Contact: jranke@uni-bremen.de
+# This file is part of the R package pfm
+
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <http://www.gnu.org/licenses/>
+
+#' Calculate a time course of relative concentrations based on an mkinmod model
+#'
+#' @import mkin
+#' @param model The degradation model to be used. Either a parent only model like
+#'   'SFO' or 'FOMC', or an mkinmod object
+#' @param DT50 The half-life. This is only used when simple exponential decline
+#'   is calculated (SFO model).
+#' @param parms The parameters used for the degradation model
+#' @param years For how many years should the degradation be predicted?
+#' @param step_days What step size in days should the output have?
+#' @param times The output times
+#' @export
+#' @author Johannes Ranke
+#' @examples
+#' head(pfm_degradation("SFO", DT50 = 10))
+pfm_degradation <- function(model = "SFO", DT50 = 1000, parms = c(k_parent_sink = log(2)/DT50),
+                            years = 1, step_days = 1,
+                            times = seq(0, years * 365, by = step_days))
+{
+  if (model %in% c("SFO", "FOMC", "DFOP", "HS", "IORE")) {
+    model <- mkinmod(parent = list(type = model))
+  }
+  initial_state = c(1, rep(0, length(model$diffs) - 1))
+  names(initial_state) <- names(model$diffs)
+  time_course <- mkinpredict(model, odeparms = parms, 
+                             odeini = initial_state,
+                             outtimes = times,
+                             solution_type = ifelse(length(model$spec) == 1,
+                                                    "analytical", "deSolve"))
+  invisible(time_course)
+}
diff --git a/R/soil_scenario_data_EFSA_2015.R b/R/soil_scenario_data_EFSA_2015.R
new file mode 100644
index 0000000..fb096bc
--- /dev/null
+++ b/R/soil_scenario_data_EFSA_2015.R
@@ -0,0 +1,40 @@
+#' Properties of the predefined scenarios from the EFSA guidance from 2015
+#' 
+#' Properties of the predefined scenarios used at Tier 1, Tier 2A and Tier 3A for the 
+#' concentration in soil as given in the EFSA guidance (2015, p. 13/14). Also, the 
+#' scenario and model adjustment factors from p. 15 and p. 17 are included.
+#' 
+#' @name soil_scenario_data_EFSA_2015
+#' @docType data
+#' @format A data frame with one row for each scenario. Row names are the scenario codes, 
+#'   e.g. CTN for the Northern scenario for the total concentration in soil. Columns are 
+#'   mostly self-explanatory. \code{rho} is the dry bulk density of the top soil.
+#' @source EFSA (European Food Safety Authority) (2015)
+#' EFSA guidance document for predicting environmental concentrations
+#' of active substances of plant protection products and transformation products of these
+#' active substances in soil. \emph{EFSA Journal} \bold{13}(4) 4093
+#' doi:10.2903/j.efsa.2015.4093
+#' @keywords datasets
+#' @examples
+#' \dontrun{
+#'   # This is the code that was used to define the data
+#'   soil_scenario_data_EFSA_2015 <- data.frame(
+#'     Zone = rep(c("North", "Central", "South"), 2),
+#'     Country = c("Estonia", "Germany", "France", "Denmark", "Czech Republik", "Spain"),
+#'     T_arit = c(4.7, 8.0, 11.0, 8.2, 9.1, 12.8),
+#'     T_arr = c(7.0, 10.1, 12.3, 9.8, 11.2, 14.7),
+#'     Texture = c("Coarse", "Coarse", "Medium fine", "Medium", "Medium", "Medium"),
+#'     f_om = c(0.118, 0.086, 0.048, 0.023, 0.018, 0.011),
+#'     theta_fc = c(0.244, 0.244, 0.385, 0.347, 0.347, 0.347),
+#'     rho = c(0.95, 1.05, 1.22, 1.39, 1.43, 1.51),
+#'     f_sce = c(3, 2, 2, 2, 1.5, 1.5),
+#'     f_mod = c(2, 2, 2, 4, 4, 4),
+#'     stringsAsFactors = FALSE,
+#'     row.names = c("CTN", "CTC", "CTS", "CLN", "CLC", "CLS")
+#'   )
+#'   save(soil_scenario_data_EFSA_2015, file = '../data/soil_scenario_data_EFSA_2015.RData')
+#' }
+#'
+#' # And this is the resulting dataframe
+#' soil_scenario_data_EFSA_2015
+NULL