Single line of generated Step12 input file partially validated

9 files changed, 303 insertions, 122 deletions

diff --git a/ChangeLog b/ChangeLog
index f749300..a64a55c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,27 @@
+commit c9bcd8e68db61515080ff377c6a04fa807337258
+Author: Johannes Ranke <jranke@uni-bremen.de>
+Date:   2017-06-17 16:36:13 +0200
+
+    Start with the generation of an input file
+
+commit e6f968cf97ed6ca9268e6098d86ba63ff2c6d2b0
+Author: Johannes Ranke <jranke@uni-bremen.de>
+Date:   2017-05-24 16:07:35 +0200
+
+    Fix for the sawtooth function for repeated applications
+    
+    For n > 2, the second application was made at 2 * i instead of i.
+
+commit 4835e20d1d08203657ab616600286ad9dfd71344
+Author: Johannes Ranke <jranke@uni-bremen.de>
+Date:   2017-05-24 15:46:44 +0200
+
+    Re-enable PELMO examples and tests
+    
+    - Add .gitattributes to make sure CRLF line endings are kept for PELMO
+    .psm files
+    - Update static docs
+
 commit 539ea37b45ddc41b36dd199f06ffe5936ab13f21
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Date:   2017-05-17 12:23:38 +0200
diff --git a/DESCRIPTION b/DESCRIPTION
index fcefaa4..6173a49 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,8 +2,8 @@ Package: pfm
 Type: Package
 Title: Utilities for Pesticide Fate Modelling
 Version: 0.4-3
-Date: 2017-05-24
-Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de", 
+Date: 2017-06-19
+Authors@R: person("Johannes Ranke", email = "jranke@uni-bremen.de",
                   role = c("aut", "cre", "cph"))
 Description: Utilities for simple calculations of predicted environmental
     concentrations (PECs) and for dealing with data from some FOCUS pesticide
diff --git a/R/PEC_sw_focus.R b/R/PEC_sw_focus.R
index e6f2689..77839ab 100644
--- a/R/PEC_sw_focus.R
+++ b/R/PEC_sw_focus.R
@@ -1,10 +1,12 @@
-#' Calculate FOCUS Step 1 PEC surface water
+#' Calculate PEC surface water at FOCUS Step 1
 #'
-#' 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
+#' This is reimplementation of the FOCUS Step 1 and 2 calculator version 3.2,
+#' authored by Michael Klein, in R. 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.
+#' this implementation. Only Step 1 PECs are calculated. However,
+#' input files are generated that are suitable as input also for Step 2
+#' to be used with the FOCUS calculator.
 #'
 #' @importFrom utils read.table
 #' @references FOCUS (2014) Generic guidance for Surface Water Scenarios (version 1.4).
@@ -40,6 +42,12 @@
 #'   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}}
+#' @param region 'n' for Northern Europe or 's' for Southern Europe. If NA, only
+#'   Step 1 PECsw are calculated
+#' @param season 'of' for October to February, 'mm' for March to May, and 'js'
+#'   for June to September. If NA, only step 1 PECsw are calculated
+#' @param interception One of 'no interception' (default), 'minimal crop cover',
+#'   'average crop cover' or 'full canopy'
 #' @param txt_file the name, and potentially the full path to the
 #'   Steps.12 input text file to which the specification of the run(s)
 #'   should be written
@@ -48,11 +56,11 @@
 #' @param append Should the input text file be appended?
 #' @examples
 #' # Parent only
-#' dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
+#' dummy_1 <- chent_focus_sw("Dummy 1", 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)
+#' new_dummy <- chent_focus_sw("New Dummy", 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,
@@ -60,6 +68,9 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
   met = NULL,
   f_drift = NA, f_rd = 0.1,
   scenario = FOCUS_Step_12_scenarios$names,
+  region = c("n", "s"),
+  season = c(NA, "of", "mm", "js"),
+  interception = "no interception",
   txt_file = "pesticide.txt", overwrite = FALSE, append = TRUE)
 {
   if (n > 1 & is.na(i)) stop("Please specify the interval i if n > 1")
@@ -67,17 +78,17 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
   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
-    # }
+    # For Step 2 we would select the reduced percentiles for multiple apps:
+    if (n <= 8) {
+      f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, as.character(n)] / 100
+    } else {
+      f_drift_2 = FOCUS_Step_12_scenarios$drift[scenario, ">8"] / 100
+    }
   }
 
+  interception = match.arg(interception)
+
   # Write to txt file if requested
-  add_line <- function(x) cat(paste0(x, "\r\n"), file = txt, append = TRUE)
-  add <- function(x) cat(paste(x, "\t"), file = txt, append = TRUE)
   if (file.exists(txt_file)) {
     if (append) {
       txt <- file(txt_file, "a")
@@ -89,52 +100,111 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
              "appending or overwriting it")
       }
     }
+  } else {
+    txt <- file(txt_file, "w")
   }
   on.exit(close(txt))
+  add_line <- function(x) cat(paste0(x, "\r\n"), file = txt, append = TRUE)
 
   # Write header to txt file
-  header <- c("Active Substance", "Compound", "Comment", "Mol mass a.i.",
-              "Mol mass met.", "Water solubility", "KOC assessed compound",
-              "KOC parent compound", "DT50", "Max. in Water",
+  header <- c("Active Substance", "Compound", "Comment",
+              "Mol mass a.i.", "Mol mass met.",
+              "Water solubility",
+              "KOC assessed compound", "KOC parent compound",
+              "DT50",
+              "Max. in Water",
               "Max. in Soil asessed compound", # we reproduce the typo...
-              "App. Rate", "Number of App.", "Time between app.",
-              "App. Type", "DT50 soil parent compound", "DT50 soil",
-              "DT50 water", "DT50 sediment", "Region / Season",
+              "App. Rate", "Number of App.", "Time between app.", "App. Type",
+              "DT50 soil parent compound", "DT50 soil",
+              "DT50 water", "DT50 sediment",
+              "Region / Season",
               "Interception class")
   add_line(paste(header, collapse = "\t"))
 
-
   if (is.null(met)) {
     compound = parent$name
     cwsat = parent$cwsat
     mw_ratio = 1
     max_soil = 1
     max_ws = 1
-    Koc = parent$Koc
+    Koc_assessed = parent$Koc
     DT50_ws = parent$DT50_ws
+    DT50_soil = parent$DT50_soil
   } else {
     compound = met$name
     cwsat = met$cwsat
     mw_ratio = met$mw / parent$mw
     max_soil = met$max_soil
     max_ws = met$max_ws
-    Koc = met$Koc
+    Koc_assessed = met$Koc
     DT50_ws = met$DT50_ws
+    DT50_soil = met$DT50_soil
   }
 
-  add(parent$name)
-  add(compound)
-  add(comment)
-  if (is.na(parent$mw)) add("-99.00")
-  else add(sprintf("%.2f", parent$mw))
-  if (is.na(met$mw)) add("-99.00")
-  else add(sprintf("%.2f", met$mw))
-  add(sprintf("%.2f", cwsat))
-  add(sprintf("%.2f", Koc))
-  if (is.null(met)) add("0.00E+00")
-  else add(sprintf("%.2f", parent$Koc))
-
-  # Rates for a single application
+  # In the Step 12 input file, the application type is the index of the drift
+  # scenarios, with numbering starting at 0
+  app_type = which(scenario == rownames(FOCUS_Step_12_scenarios$drift)) - 1
+  # The interception class is a number from 1.00 (no interception) to 4.00 (full canopy)
+  int_class = which(interception == colnames(FOCUS_Step_12_scenarios$interception))
+
+  # We calculate Step 2 if a season is specified
+  # Interception and region/season code
+  if (is.na(season[1])) {
+    f_int = 0
+    reg_sea = 0
+  } else {
+    f_int = FOCUS_Step_12_scenarios$interception[scenario, interception]
+
+    # the code for region / season is the sum of the region code
+    reg_code = switch(region,
+      "n" = 0,
+      "s" = 3)
+    # and the season code
+    sea_code = switch(season,
+      "of" = 0,
+      "mm" = 1,
+      "js" = 2)
+    reg_sea = reg_code + sea_code
+  }
+
+  run_txt <- c(
+    ai = parent$name, compound = compound,
+    comment = comment)
+
+  run_numeric <- c(
+    mw_ai = NA, mw_met = NA,
+    cwsat = cwsat, Koc_assessed = Koc_assessed, Koc_parent = 0,
+    DT50_ws = DT50_ws,
+    max_ws = 0, max_soil = 0,
+    rate = rate, n = n, i = if(is.na(i)) 0 else i,
+    app_type = app_type,
+    DT50_soil_parent = 0, DT50_soil = 0, DT50_water = 0, DT50_sediment = 0,
+    reg_sea = reg_sea, int_class = int_class)
+
+  if (is.null(met)) {
+    run_numeric["DT50_soil"] = parent$DT50_soil
+    run_numeric["DT50_water"] = parent$DT50_water
+    run_numeric["DT50_sediment"] = parent$DT50_sediment
+  } else {
+    run_numeric["mw_ai"] = parent$mw
+    run_numeric["mw_met"] = met$mw
+    run_numeric["Koc_parent"] = parent$Koc
+    run_numeric["DT50_soil"] = met$DT50_soil
+    run_numeric["DT50_soil_parent"] = parent$DT50_soil
+    run_numeric["DT50_water"] = met$DT50_water
+    run_numeric["DT50_sediment"] = met$DT50_sediment
+  }
+
+  print_numeric <- function(x) {
+    ifelse(is.na(x), "  -99.00",
+      ifelse(x < 0.01,
+        sprintf("%8.2E", x),
+        sprintf("%8.2f", x)))
+  }
+  run_line <- paste(c(run_txt, print_numeric(run_numeric)), collapse = "\t")
+  add_line(run_line)
+
+  # Rates for a single application (suffix _s)
   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
@@ -156,6 +226,7 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
   input_rd_s = f_rd * eq_rate_rd_s_tot * ratio_field_wb / 10
   input_rd = n * input_rd_s
 
+  # Step 1 accumulation in the aquatic system
   # No accumulation between multiple applications if 3 * DT50 < i
   if (n > 1 && (3 * DT50_ws < i)) {
     input_drift = input_drift_s
@@ -165,13 +236,19 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
     input_rd = n * input_rd_s
   }
 
+  # Step 2 accumulation in soil
+  # Use interception and degradation in soil for Step 2
+  k_soil = log(2)/DT50_soil
+  f_deg_soil = exp(- 4 * k_soil) * # Four days of degradation
+    (1 - exp(-n * i * k_soil)) / (1 - exp(-i * k_soil)) # multiple applications
+
   # Fraction of compound entering the water phase via runoff/drainage
   depth_sw  = 0.3  # m
   depth_sed = 0.05 # m
   depth_sed_eff = 0.01 # m, only relevant for sorption
   rho_sed = 0.8 # kg/L
   f_OC = 0.05 # 5% organic carbon in sediment
-  f_rd_sw = depth_sw / (depth_sw + (depth_sed_eff * rho_sed * f_OC * Koc))
+  f_rd_sw = depth_sw / (depth_sw + (depth_sed_eff * rho_sed * f_OC * Koc_assessed))
   f_rd_sed = 1 - f_rd_sw
 
   # Initial PECs
@@ -237,6 +314,9 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
 #' @param name Length one character vector containing the name
 #' @param cwsat Water solubility in mg/L
 #' @param DT50_ws Half-life in water/sediment systems in days
+#' @param DT50_soil Half-life in soil in days
+#' @param DT50_water Half-life in water in days (Step 2)
+#' @param DT50_sediment Half-life in sediment in days (Step 2)
 #' @param Koc Partition coefficient between organic carbon and water
 #'   in L/kg.
 #' @param mw Molar weight in g/mol.
@@ -244,8 +324,11 @@ PEC_sw_focus <- function(parent, rate, n = 1, i = NA,
 #' @param max_ws Maximum observed fraction (dimensionless) in water/sediment
 #'   systems
 #' @return A list with the substance specific properties
-chent_focus_sw <- function(name, Koc, DT50_ws = NA, cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1)
+chent_focus_sw <- function(name, Koc, DT50_ws = NA, DT50_soil = NA,
+                           DT50_water = NA, DT50_sediment = NA,
+                           cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1)
 {
-  list(Koc = Koc, DT50_ws = DT50_ws, cwsat = cwsat,
+  list(name = name, Koc = Koc, DT50_ws = DT50_ws, DT50_soil = DT50_soil,
+       DT50_water = DT50_water, DT50_sediment = DT50_sediment, cwsat = cwsat,
        mw = mw, max_soil = max_soil, max_ws = max_ws)
 }
diff --git a/build.log b/build.log
index 4199788..2b7c5e6 100644
--- a/build.log
+++ b/build.log
@@ -1,6 +1,9 @@
 * checking for file ‘./DESCRIPTION’ ... OK
 * preparing ‘pfm’:
 * checking DESCRIPTION meta-information ... OK
+* excluding invalid files
+Subdirectory 'R' contains invalid file names:
+  ‘pesticide.txt’ ‘xxx.txt’
 * checking for LF line-endings in source and make files
 * checking for empty or unneeded directories
 * looking to see if a ‘data/datalist’ file should be added
diff --git a/inst/testdata/Steps_12_pesticide.txt b/inst/testdata/Steps_12_pesticide.txt
new file mode 100644
index 0000000..a2c3768
--- /dev/null
+++ b/inst/testdata/Steps_12_pesticide.txt
@@ -0,0 +1,12 @@
+Active Substance	Compound	Comment	Mol mass a.i.	Mol mass met.	Water solubility	KOC assessed compound	KOC parent compound	DT50	Max. in Water	Max. in Soil asessed compound	App. Rate	Number of App.	Time between app.	App. Type	DT50 soil parent compound	DT50 soil	DT50 water	DT50 sediment	Region / Season	Interception class
+Dummy 1	Dummy 1	Potatoes, Southern Europe, spring, 1 app/season, soil incorporation	  -99.00	  -99.00	 6000.00	  344.80	0.00E+00	    6.00	  100.00	  100.00	 3000.00	    1.00	0.00E+00	   28.00	 1000.00	    6.00	    6.00	    6.00	    4.00	    1.00
+Dummy 2	Dummy 2	Maize, Southern Europe, spring, 1 app/season	  100.00	  -99.00	   30.00	  110.00	0.00E+00	   26.00	0.00E+00	0.00E+00	 1000.00	    1.00	0.00E+00	    8.00	 1000.00	   56.00	   26.00	   26.00	    4.00	    1.00
+Dummy 3	Dummy 3	Winter wheat,, Northern  Europe, spring, 1 app/season	0.00E+00	  -99.00	  620.00	    1.00	0.00E+00	    1.50	0.00E+00	0.00E+00	 1000.00	    1.00	0.00E+00	    1.00	 1000.00	    4.00	    1.50	    1.50	    1.00	    1.00
+Dummy 4	Dummy 4	Apples, Southern  Europe, spring, 3 app./season, 14 d int, orchards	0.00E+00	  -99.00	2.00E-03	  970.00	0.00E+00	    4.00	0.00E+00	0.00E+00	    7.50	    3.00	   14.00	   12.00	 1000.00	   19.00	    4.00	    4.00	    4.00	    1.00
+Dummy 5	Dummy 5	Vines, Northern Europe, spring, 5 app/seaon 14 d int.	0.00E+00	  -99.00	    1.15	  860.00	0.00E+00	  118.00	0.00E+00	0.00E+00	   75.00	    5.00	   14.00	   23.00	 1000.00	  250.00	    6.00	  118.00	    1.00	    1.00
+Dummy 6	Dummy 6	Cereals, Northern Europe, spring, 1 app/season	0.00E+00	  -99.00	 6000.00	   66.00	0.00E+00	   24.00	0.00E+00	0.00E+00	  400.00	    1.00	0.00E+00	0.00E+00	 1000.00	   28.00	   24.00	   24.00	    1.00	    1.00
+Dummy 6	Dummy 6	Cereals, Southern Europe, spring, 1 app/seaon	0.00E+00	  -99.00	  200.00	   66.00	0.00E+00	   24.00	0.00E+00	0.00E+00	  200.00	    1.00	0.00E+00	0.00E+00	 1000.00	   28.00	   24.00	   24.00	    4.00	    1.00
+Dummy 7	Dummy 7	Vines, Southern Europe, spring, 4 app/seaon 14 d int.	0.00E+00	  -99.00	    2.60	  500.00	0.00E+00	   28.00	0.00E+00	0.00E+00	  750.00	    4.00	   14.00	   23.00	 1000.00	   50.00	    2.50	   28.00	    4.00	    1.00
+New Dummy	Metabolite M1	Soil Metabolite	  250.00	  100.00	  100.00	   50.00	  100.00	  100.00	0.00E+00	   50.00	 1000.00	    1.00	0.00E+00	    1.00	   10.00	   20.00	   10.00	  100.00	0.00E+00	    1.00
+New Dummy	Metabolite M2	Water Metabolite	  250.00	  100.00	  100.00	   50.00	  100.00	  100.00	   50.00	0.00E+00	 1000.00	    1.00	0.00E+00	    1.00	   10.00	   20.00	   10.00	  100.00	0.00E+00	    1.00
+New Dummy	Metabolite M3	Soil and Water Metabolite	  250.00	  100.00	  100.00	   50.00	  100.00	  100.00	   50.00	   50.00	 1000.00	    1.00	0.00E+00	    1.00	   10.00	   20.00	   10.00	  100.00	0.00E+00	    1.00
diff --git a/man/PEC_sw_focus.Rd b/man/PEC_sw_focus.Rd
index a5628d7..e9f638d 100644
--- a/man/PEC_sw_focus.Rd
+++ b/man/PEC_sw_focus.Rd
@@ -2,11 +2,13 @@
 % Please edit documentation in R/PEC_sw_focus.R
 \name{PEC_sw_focus}
 \alias{PEC_sw_focus}
-\title{Calculate FOCUS Step 1 PEC surface water}
+\title{Calculate PEC surface water at FOCUS Step 1}
 \usage{
-PEC_sw_focus(parent, rate, n = 1, i = NA, met = NULL, f_drift = NA,
-  f_rd = 0.1, scenario = FOCUS_Step_12_scenarios$names,
-  txt_file = "pesticide.txt", overwrite = FALSE, append = TRUE)
+PEC_sw_focus(parent, rate, n = 1, i = NA, comment = "", met = NULL,
+  f_drift = NA, f_rd = 0.1, scenario = FOCUS_Step_12_scenarios$names,
+  region = c("n", "s"), season = c(NA, "of", "mm", "js"),
+  interception = "no interception", txt_file = "pesticide.txt",
+  overwrite = FALSE, append = TRUE)
 }
 \arguments{
 \item{parent}{A list containing substance specific parameters, e.g.
@@ -19,6 +21,8 @@ applications are given explicitly}
 
 \item{i}{The application interval}
 
+\item{comment}{A comment for the input file}
+
 \item{met}{A list containing metabolite specific parameters.  e.g.
 conveniently generated by \code{\link{chent_focus_sw}}.  If not NULL,
 the PEC is calculated for this compound, not the parent.}
@@ -35,6 +39,15 @@ parent or a metabolite}
 \item{scenario}{The name of the scenario. Must be one of the scenario
 names given in \code{\link{FOCUS_Step_12_scenarios}}}
 
+\item{region}{'n' for Northern Europe or 's' for Southern Europe. If NA, only
+Step 1 PECsw are calculated}
+
+\item{season}{'of' for October to February, 'mm' for March to May, and 'js'
+for June to September. If NA, only step 1 PECsw are calculated}
+
+\item{interception}{One of 'no interception' (default), 'minimal crop cover',
+'average crop cover' or 'full canopy'}
+
 \item{txt_file}{the name, and potentially the full path to the
 Steps.12 input text file to which the specification of the run(s)
 should be written}
@@ -45,11 +58,13 @@ should be written}
 \item{append}{Should the input text file be appended?}
 }
 \description{
-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
+This is reimplementation of the FOCUS Step 1 and 2 calculator version 3.2,
+authored by Michael Klein, in R. 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.
+this implementation. Only Step 1 PECs are calculated. However,
+input files are generated that are suitable as input also for Step 2
+to be used with the FOCUS calculator.
 }
 \note{
 The formulas for input to the waterbody via runoff/drainage of the
@@ -63,11 +78,11 @@ Step 2 is not implemented.
 }
 \examples{
 # Parent only
-dummy_1 <- chent_focus_sw(cwsat = 6000, DT50_ws = 6, Koc = 344.8)
+dummy_1 <- chent_focus_sw("Dummy 1", 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)
+new_dummy <- chent_focus_sw("New Dummy", 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)
 }
diff --git a/man/chent_focus_sw.Rd b/man/chent_focus_sw.Rd
index 945c343..8df01ab 100644
--- a/man/chent_focus_sw.Rd
+++ b/man/chent_focus_sw.Rd
@@ -4,18 +4,26 @@
 \alias{chent_focus_sw}
 \title{Create a chemical compound object for FOCUS Step 1 calculations}
 \usage{
-chent_focus_sw(Koc, DT50_ws = NA, cwsat = 1000, mw = NA, max_soil = 1,
-  max_ws = 1)
+chent_focus_sw(name, Koc, DT50_ws = NA, DT50_soil = NA, DT50_water = NA,
+  DT50_sediment = NA, cwsat = 1000, mw = NA, max_soil = 1, max_ws = 1)
 }
 \arguments{
+\item{name}{Length one character vector containing the name}
+
 \item{Koc}{Partition coefficient between organic carbon and water
 in L/kg.}
 
 \item{DT50_ws}{Half-life in water/sediment systems in days}
 
+\item{DT50_soil}{Half-life in soil in days}
+
+\item{DT50_water}{Half-life in water in days (Step 2)}
+
+\item{DT50_sediment}{Half-life in sediment in days (Step 2)}
+
 \item{cwsat}{Water solubility in mg/L}
 
-\item{mw}{Molar weight in g/mol}
+\item{mw}{Molar weight in g/mol.}
 
 \item{max_soil}{Maximum observed fraction (dimensionless) in soil}
 
diff --git a/test.log b/test.log
index 2963aae..af3bddb 100644
--- a/test.log
+++ b/test.log
@@ -15,7 +15,7 @@ Simple PEC soil calculations: ...........
 Simple PEC surface water calculations with drift entry: ..
 Create PELMO runs from psm files and execute them: ...........................................................................................................................................................................................................................
 Actual and time weighted average concentrations for SFO kinetics: .
-FOCUS Step 1 calculations: .........
+FOCUS Step 1 and 2 calculations: ...........
 Read and analyse TOXSWA cwa files: .......
 UK drainage PEC calculations: ............
 
diff --git a/tests/testthat/test_step_1.R b/tests/testthat/test_step_1.R
index 0c8710b..2a65f36 100644
--- a/tests/testthat/test_step_1.R
+++ b/tests/testthat/test_step_1.R
@@ -1,125 +1,161 @@
-context("FOCUS Step 1 calculations")
+context("FOCUS Step 1 and 2 calculations")
 
-t_out <- c(0, 1, 2, 4) # Checking the first four days is sufficient for Step 1
 test_txt <- readLines(
   system.file("testdata/Steps_12_pesticide.txt", package = "pfm")
 )
 
+# Define test compounds as in pesticide.txt
+dummy_1 <- chent_focus_sw("Dummy 1", cwsat = 6000, DT50_ws = 6, DT50_soil = 6, Koc = 344.8,
+                          DT50_water = 6, DT50_sediment = 6)
+dummy_2 <- chent_focus_sw("Dummy 2", cwsat = 30, DT50_ws = 26, Koc = 110)
+dummy_4 <- chent_focus_sw("Dummy 4", cwsat = 2e-3, DT50_ws = 4, Koc = 970)
+dummy_5 <- chent_focus_sw("Dummy 5", cwsat = 1.15, DT50_ws = 118, Koc = 860)
+dummy_7 <- chent_focus_sw("Dummy 7", cwsat = 2.60, DT50_ws = 28, Koc = 500)
+new_dummy <- chent_focus_sw("New Dummy", mw = 250, Koc = 100)
+M1 <- chent_focus_sw("M1", mw = 100, cwsat = 100, DT50_ws = 100,
+                     Koc = 50, max_ws = 0, max_soil = 0.5)
+M2 <- chent_focus_sw("M2", mw = 100, cwsat = 100, DT50_ws = 100,
+                     Koc = 50, max_ws = 0.5, max_soil = 0)
+
+# When we compare the generated input file with the test file,
+# we can ignore some fields
+field_index <- c(ai = 1, compound = 2, comment = 3,
+  mw_ai = 4, mw_met = 5,
+  cwsat = 6, Koc_assessed = 7,
+  Koc_parent = 8,
+  DT50_ws = 9,
+  max_ws = 10, max_soil = 11,
+  rate = 12, n = 13, i = 14, app_type = 15,
+  DT50_soil_parent = 16, DT50_soil = 17, DT50_water = 18, DT50_sediment = 19,
+  reg_sea = 20, int_class = 21)
+field_index_parent <- field_index[-c(4:5, 8, 10:11, 16)]
+
+t_out_1 <- c(0, 1, 2, 4) # Checking the first four days is sufficient for Step 1
+PEC_template_1 <- matrix(NA, nrow = length(t_out_1), ncol = 4,
+    dimnames = list(Time = t_out_1, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+
+t_out_2 <- c(0, 1, 2, 4, 7, 14, 21, 28, 42, 50, 100) # We read in text from rtf reports for Step 2
+
 test_that("Results of Steps 1/2 calculator for Dummy 1 are reproduced", {
-  dummy_1 <- chent_focus_sw("Dummy 1", cwsat = 6000, DT50_ws = 6, Koc = 344.8)
-  res_dummy_1 <- PEC_sw_focus(dummy_1, 3000,
+  res_step_1_1 <- PEC_sw_focus(dummy_1, 3000,
     comment = "Potatoes, Southern Europe, spring, 1 app/season, soil incorporation",
-    f_drift = 0,
+    scenario = "no drift (incorp or seed trtmt)",
+    region = "s", season = "mm",
     append = FALSE, overwrite = TRUE)
 
   pest_txt <- readLines("pesticide.txt")
-  expect_equal(test_txt[1], pest_txt[1])
-  strsplit(test_txt[2], "\t")[[1]]
-
-  PEC_orig_1 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
-
-  PEC_orig_1[, "PECsw"]  = c(685.06, 610.32, 543.73, 431.56)
-  PEC_orig_1[, "TWAECsw"]  = c(NA, 647.69, 612.03, 548.76)
-  PEC_orig_1[, "PECsed"]  = c(2.36, 2.1, 1.87, 1.49) * 1e3
-  PEC_orig_1[, "TWAECsed"]  = c(NA, 2.23e3, 2.11e3, 1.89e3)
-
-  expect_equal(res_dummy_1$PEC[1:4, c(1, 2)], PEC_orig_1[, c(1, 2)], tolerance = 0.01, scale = 1)
-  expect_equal(res_dummy_1$PEC[1:4, c(3, 4)], PEC_orig_1[, c(3, 4)], tolerance = 10, scale = 1)
+  expect_equal(test_txt[1], pest_txt[1]) # Header
+  test_1 <- strsplit(test_txt[2], "\t")[[1]][field_index_parent]
+  pest_1 <- strsplit(pest_txt[2], "\t")[[1]][field_index_parent]
+  expect_equal(test_1, pest_1) # Parent fields
+
+  PEC_step_1_1 <- PEC_template_1
+  PEC_step_1_1[, "PECsw"]  = c(685.06, 610.32, 543.73, 431.56)
+  PEC_step_1_1[, "TWAECsw"]  = c(NA, 647.69, 612.03, 548.76)
+  PEC_step_1_1[, "PECsed"]  = c(2.36, 2.1, 1.87, 1.49) * 1e3
+  PEC_step_1_1[, "TWAECsed"]  = c(NA, 2.23e3, 2.11e3, 1.89e3)
+
+  expect_equal(res_step_1_1$PEC[1:4, c(1, 2)], PEC_step_1_1[, c(1, 2)], tolerance = 0.01, scale = 1)
+  expect_equal(res_step_1_1$PEC[1:4, c(3, 4)], PEC_step_1_1[, c(3, 4)], tolerance = 10, scale = 1)
+
+  # This is pasted from the file "Dummy 1 step 2.rtf" generated with Steps12 version 3.2 from 15/05/2017
+  PEC_step_2_1_raw <- read.table(text = "0 172.6235  NA 595.2057  NA
+1 153.7900  163.2067  530.2680  562.7368
+2 137.0113  154.3037  472.4151  532.0392
+4 108.7460  138.3873  374.9561  477.1595
+7  76.8950  118.5090  265.1340  408.6191
+14   34.2528   85.6494  118.1038  295.3191
+21   15.2579   64.9380   52.6092  223.9062
+28    6.7966   51.3222   23.4348  176.9589
+42    1.3486   35.3389    4.6500  121.8484
+50    0.5352   29.8256    1.8454  102.8388
+100   0.0017   14.9591    0.0057   51.5788")
+  PEC_step_2_1 = PEC_step_2_1_raw[, 2:5]
+  dimnames(PEC_step_2_1) = list(Time = t_out_2,
+    type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed"))
+
+  # Step 2 is not implemented.
 })
 
 test_that("Results of Steps 1/2 calculator for Dummy 2 are reproduced", {
-  dummy_2 <- chent_focus_sw(cwsat = 30, DT50_ws = 26, Koc = 110)
   res_dummy_2 <- PEC_sw_focus(dummy_2, 1000)
 
-  PEC_orig_2 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_2 = PEC_template_1
 
-  PEC_orig_2[, "PECsw"]  = c(299.89, 290.86, 283.21, 268.50)
-  PEC_orig_2[, "TWAECsw"]  = c(NA, 295.38, 291.20, 283.49)
-  PEC_orig_2[, "PECsed"]  = c(319.77, 319.95, 311.53, 295.35)
-  PEC_orig_2[, "TWAECsed"]  = c(NA, 319.86, 317.79, 310.58)
+  PEC_step_1_2[, "PECsw"]  = c(299.89, 290.86, 283.21, 268.50)
+  PEC_step_1_2[, "TWAECsw"]  = c(NA, 295.38, 291.20, 283.49)
+  PEC_step_1_2[, "PECsed"]  = c(319.77, 319.95, 311.53, 295.35)
+  PEC_step_1_2[, "TWAECsed"]  = c(NA, 319.86, 317.79, 310.58)
 
-  expect_equal(res_dummy_2$PEC[1:4, ], PEC_orig_2[, ], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_2$PEC[1:4, ], PEC_step_1_2[, ], tolerance = 0.01, scale = 1)
 })
 
 test_that("Results of Steps 1/2 calculator for Dummy 4 are reproduced", {
-  dummy_4 <- chent_focus_sw(cwsat = 2e-3, DT50_ws = 4, Koc = 970)
   res_dummy_4 <- PEC_sw_focus(dummy_4, 7.5, n = 3, i = 14,
                               scenario = "pome / stone fruit, early")
 
-  PEC_orig_4 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_4 = PEC_template_1
 
-  PEC_orig_4[, "PECsw"]  = c(1.82, 1.18, 1.00, 0.70)
-  PEC_orig_4[, "TWAECsw"]  = c(NA, 1.50, 1.29, 1.07)
-  PEC_orig_4[, "PECsed"]  = c(10.57, 11.49, 9.66, 6.83)
-  PEC_orig_4[, "TWAECsed"]  = c(NA, 11.03, 10.79, 9.48)
+  PEC_step_1_4[, "PECsw"]  = c(1.82, 1.18, 1.00, 0.70)
+  PEC_step_1_4[, "TWAECsw"]  = c(NA, 1.50, 1.29, 1.07)
+  PEC_step_1_4[, "PECsed"]  = c(10.57, 11.49, 9.66, 6.83)
+  PEC_step_1_4[, "TWAECsed"]  = c(NA, 11.03, 10.79, 9.48)
 
-  expect_equal(res_dummy_4$PEC[1:4, ], PEC_orig_4[, ], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_4$PEC[1:4, ], PEC_step_1_4[, ], tolerance = 0.01, scale = 1)
 })
 
 test_that("Results of Steps 1/2 calculator for Dummy 5 are reproduced", {
-  dummy_5 <- chent_focus_sw(cwsat = 1.15, DT50_ws = 118, Koc = 860)
   res_dummy_5 <- PEC_sw_focus(dummy_5, 75, n = 5, i = 14,
                               scenario = "vines, early")
 
-  PEC_orig_5 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_5 = PEC_template_1
 
-  PEC_orig_5[, "PECsw"]  = c(61.60, 59.45, 59.10, 58.41)
-  PEC_orig_5[, "TWAECsw"]  = c(NA, 60.53, 59.90, 59.33)
-  PEC_orig_5[, "PECsed"]  = c(500.78, 511.28, 508.29, 502.35)
-  PEC_orig_5[, "TWAECsed"]  = c(NA, 506.03, 507.90, 506.61)
+  PEC_step_1_5[, "PECsw"]  = c(61.60, 59.45, 59.10, 58.41)
+  PEC_step_1_5[, "TWAECsw"]  = c(NA, 60.53, 59.90, 59.33)
+  PEC_step_1_5[, "PECsed"]  = c(500.78, 511.28, 508.29, 502.35)
+  PEC_step_1_5[, "TWAECsed"]  = c(NA, 506.03, 507.90, 506.61)
 
-  expect_equal(res_dummy_5$PEC[1:4, ], PEC_orig_5[, ], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_5$PEC[1:4, ], PEC_step_1_5[, ], tolerance = 0.01, scale = 1)
 })
 
 test_that("Results of Steps 1/2 calculator for Dummy 7 are reproduced", {
-  dummy_7 <- chent_focus_sw(cwsat = 2.60, DT50_ws = 28, Koc = 500)
   res_dummy_7 <- PEC_sw_focus(dummy_7, 750, n = 4, i = 14,
                               scenario = "vines, early")
 
-  PEC_orig_7 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_7 = PEC_template_1
 
-  PEC_orig_7[, "PECsw"]  = c(626.99, 601.13, 586.43, 558.10)
-  PEC_orig_7[, "TWAECsw"]  = c(NA, 614.06, 603.90, 588.03)
-  PEC_orig_7[, "PECsed"]  = c(3.0, 3.01, 2.93, 2.79) * 1e3
-  PEC_orig_7[, "TWAECsed"]  = c(NA, 3.01e3, 2.99e3, 2.92e3)
+  PEC_step_1_7[, "PECsw"]  = c(626.99, 601.13, 586.43, 558.10)
+  PEC_step_1_7[, "TWAECsw"]  = c(NA, 614.06, 603.90, 588.03)
+  PEC_step_1_7[, "PECsed"]  = c(3.0, 3.01, 2.93, 2.79) * 1e3
+  PEC_step_1_7[, "TWAECsed"]  = c(NA, 3.01e3, 2.99e3, 2.92e3)
 
-  expect_equal(res_dummy_7$PEC[1:4, c(1, 2)], PEC_orig_7[, c(1, 2)], tolerance = 0.01, scale = 1)
-  expect_equal(res_dummy_7$PEC[1:4, c(3, 4)], PEC_orig_7[, c(3, 4)], tolerance = 10, scale = 1)
+  expect_equal(res_dummy_7$PEC[1:4, c(1, 2)], PEC_step_1_7[, c(1, 2)], tolerance = 0.01, scale = 1)
+  expect_equal(res_dummy_7$PEC[1:4, c(3, 4)], PEC_step_1_7[, c(3, 4)], tolerance = 10, scale = 1)
 })
 
 test_that("Results of Steps 1/2 calculator for New Dummy (M1-M3) are reproduced", {
-  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)
   res_M1 <- PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter",
                          met = M1)
 
-  PEC_orig_M1 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_M1 = PEC_template_1
 
-  PEC_orig_M1[, "PECsw"]  = c(62.5, 62.07, 61.64, 60.79)
-  PEC_orig_M1[, "TWAECsw"]  = c(NA, 62.28, 62.07, 61.64)
-  PEC_orig_M1[, "PECsed"]  = c(31.25, 31.03, 30.82, 30.40)
-  PEC_orig_M1[, "TWAECsed"]  = c(NA, 31.14, 31.03, 30.82)
+  PEC_step_1_M1[, "PECsw"]  = c(62.5, 62.07, 61.64, 60.79)
+  PEC_step_1_M1[, "TWAECsw"]  = c(NA, 62.28, 62.07, 61.64)
+  PEC_step_1_M1[, "PECsed"]  = c(31.25, 31.03, 30.82, 30.40)
+  PEC_step_1_M1[, "TWAECsed"]  = c(NA, 31.14, 31.03, 30.82)
 
-  expect_equal(res_M1$PEC[1:4, ], PEC_orig_M1[, ], tolerance = 0.01, scale = 1)
+  expect_equal(res_M1$PEC[1:4, ], PEC_step_1_M1[, ], tolerance = 0.01, scale = 1)
 
-  M2 <- chent_focus_sw(mw = 100, cwsat = 100, DT50_ws = 100, Koc = 50, max_ws = 0.5, max_soil = 0)
   res_M2 <- PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter",
                          met = M2)
 
-  PEC_orig_M2 = matrix(NA, nrow = length(t_out), ncol = 4,
-    dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))
+  PEC_step_1_M2 = PEC_template_1
 
-  PEC_orig_M2[, "PECsw"]  = c(64.34, 63.78, 63.34, 62.47)
-  PEC_orig_M2[, "TWAECsw"]  = c(NA, 64.06, 63.81, 63.36)
-  PEC_orig_M2[, "PECsed"]  = c(31.25, 31.89, 31.67, 31.23)
-  PEC_orig_M2[, "TWAECsed"]  = c(NA, 31.57, 31.68, 31.56)
+  PEC_step_1_M2[, "PECsw"]  = c(64.34, 63.78, 63.34, 62.47)
+  PEC_step_1_M2[, "TWAECsw"]  = c(NA, 64.06, 63.81, 63.36)
+  PEC_step_1_M2[, "PECsed"]  = c(31.25, 31.89, 31.67, 31.23)
+  PEC_step_1_M2[, "TWAECsed"]  = c(NA, 31.57, 31.68, 31.56)
 
-  expect_equal(res_M2$PEC[1:4, ], PEC_orig_M2[, ], tolerance = 0.01, scale = 1)
+  expect_equal(res_M2$PEC[1:4, ], PEC_step_1_M2[, ], tolerance = 0.01, scale = 1)
 })
 unlink("pesticide.txt")