#' Calculate PEC surface water at FOCUS Step 1
#'
#' This is a 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. Only Step 1 PECs are calculated. However,
#' input files can be generated that are suitable as input for 
#' the FOCUS calculator.
#'
#' @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. As one would
#'   expect, this appears to be (as we get the same results) calculated by
#'   multiplying the application rate with the molar weight
#'   correction and the formation fraction in water/sediment systems.
#' @note Step 2 is not implemented.
#' @export
#' @param parent A list containing substance specific parameters, e.g.
#'   conveniently generated by [chent_focus_sw].
#' @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 comment A comment for the input file
#' @param met A list containing metabolite specific parameters.  e.g.
#'   conveniently generated by [chent_focus_sw].  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
#'   [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 [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 met_form_water Should the metabolite formation in water be taken into
#'   account? This can be switched off to check the influence and to compare
#'   with previous versions of the Steps 12 calculator
#' @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
#' @param overwrite Should an existing file a the location specified in
#'   \code{txt_file} be overwritten? Only takes effect if append is FALSE.
#' @param append Should the input text file be appended, if it exists?
#' @examples
#' # Parent only
#' 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("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)
#' PEC_sw_focus(new_dummy, 1000, scenario = "cereals, winter", met = M1)
PEC_sw_focus <- function(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 = c("no interception", "minimal crop cover",
                   "average crop cover", "full canopy"),
  met_form_water = TRUE,
  txt_file = "pesticide.txt", overwrite = FALSE, append = FALSE)
{
  if (n > 1 & is.na(i)) stop("Please specify the interval i if n > 1")

  scenario = match.arg(scenario)

  if (is.na(f_drift)) {
    f_drift = FOCUS_Step_12_scenarios$drift[scenario, "1"] / 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)

  if (append | overwrite) {
    # Write to txt file if requested
    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",
                "Interception class")
    add_line <- function(x) {
      cat(paste0(x, "\r\n"), file = txt, append = TRUE)
    }
    if (file.exists(txt_file)) {
      if (append) {
        txt <- file(txt_file, "a")
      } else {
        if (overwrite) {
          txt <- file(txt_file, "w")
          add_line(paste(header, collapse = "\t"))
        }
      }
    } else {
      txt <- file(txt_file, "w")
      add_line(paste(header, collapse = "\t"))
    }
    on.exit(close(txt))
  }

  region = match.arg(region)
  season = match.arg(season)
  interception = match.arg(interception)

  if (is.null(met)) {
    compound = parent$name
    cwsat = parent$cwsat
    mw_ratio = 1
    max_soil = 1
    max_ws = 1
    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_assessed = met$Koc
    DT50_ws = met$DT50_ws
    DT50_soil = met$DT50_soil
  }

  # 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
  }

  scenario_safe <- sub(" /", " or", scenario)
  if (is.null(met)) {
    name_input <- paste(parent$name, scenario_safe, region, season)
  } else {
    name_input <- paste(met$name, scenario_safe, region, season)
  }
  if (comment != "") name_input = paste(name_input, comment)

  run_txt <- c(
    ai = name_input, compound = name_input,
    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["max_soil"] = 100 * met$max_soil
    run_numeric["max_ws"] = 100 * met$max_ws
    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")
  if (append | overwrite) {
    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
  # Metabolite formed in water (this part is not documented in the Help files
  # of the Steps 1/2 calculator):
  if (!is.null(met)) {
    if (met_form_water) {
      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
    }
  } 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

  # 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_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
    input_rd = input_rd_s
  } else {
    input_drift = n * input_drift_s
    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_assessed))
  f_rd_sed = 1 - f_rd_sw

  # Initial PECs
  PEC_sw_0 = (input_drift + input_rd * f_rd_sw) / depth_sw        # µg/L
  PEC_sed_0 = (input_rd * f_rd_sed) / (depth_sed * rho_sed)       # µg/kg

  # Initial PECs when assuming partitioning also of drift input
  PEC_sw_0_part = (input_drift + input_rd) * f_rd_sw / depth_sw                  # µg/L
  PEC_sed_0_part = (input_drift + input_rd) * f_rd_sed / (depth_sed * rho_sed)   # µg/kg

  t_out = c(0, 1, 2, 4, 7, 14, 21, 28, 42, 50, 100)
  PEC = matrix(NA, nrow = length(t_out), ncol = 4,
               dimnames = list(Time = t_out, type = c("PECsw", "TWAECsw", "PECsed", "TWAECsed")))

  PEC["0", "PECsw"]  = PEC_sw_0
  PEC["0", "PECsed"] = PEC_sed_0

  # Degradation after partitioning according to Koc
  k_ws = log(2)/DT50_ws
  PEC[as.character(t_out[-1]), "PECsw"]  = PEC_sw_0_part  * exp( - k_ws * t_out[-1])
  PEC[as.character(t_out[-1]), "PECsed"] = PEC_sed_0_part * exp( - k_ws * t_out[-1])

  # TWA concentrations
  PEC_sw_1 = PEC["1", "PECsw"]
  PEC_sed_1 = PEC["1", "PECsed"]
  TWAEC_sw_1 = (PEC_sw_0 + PEC_sw_1) / 2
  TWAEC_sed_1 = (PEC_sed_0 + PEC_sed_1) / 2

  TWAEC_after = function(TWAEC_1, PEC_1, t) {
    (TWAEC_1 + PEC_1 * (1 - exp(- k_ws * (t - 1))) / k_ws) / t
  }

  PEC["1", "TWAECsw"] = TWAEC_sw_1
  PEC["1", "TWAECsed"] = TWAEC_sed_1

  PEC[as.character(t_out[-1]), "TWAECsw"] = TWAEC_after(TWAEC_sw_1, PEC_sw_1, t_out[-1])
  PEC[as.character(t_out[-1]), "TWAECsed"] = TWAEC_after(TWAEC_sed_1, PEC_sed_1, t_out[-1])

  # Check if PEC_sw_max is above water solubility
  PEC_sw_max = max(PEC[, "PECsw"])
  if (!is.na(PEC_sw_max) & 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 a chemical compound object for FOCUS Step 1 calculations
#'
#' @export
#' @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.
#' @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(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(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)
}