Separate out PELMO utilities into rPELMO

2 files changed, 0 insertions, 571 deletions

diff --git a/R/FOCUS_PELMO_data.R b/R/FOCUS_PELMO_data.R
deleted file mode 100644
index 89fe5f6..0000000
--- a/R/FOCUS_PELMO_data.R
+++ /dev/null
@@ -1,90 +0,0 @@
-#' FOCUS PELMO crop names
-#'
-#' A named character vector with the crop names used in the PELMO 5.5.3 GUI.
-#' For the names, three letter codes were constructed by generally taking the
-#' first three letters in lower case. Only when there is an expression in
-#' parentheses, the first letter of this expression becomes the first letter
-#' in the three letter code, i.e. 'Peas (animals)' has the code \code{ape}.
-#'
-#' @docType data
-#' @export
-#' @examples
-#' print(FOCUS_PELMO_crops)
-FOCUS_PELMO_crops = c(
-  app = "Apples",
-  gra = "Grass and alfalfa",
-  pot = "Potatoes",
-  sug = "Sugar beets",
-  win = "Winter cereals",
-  fbe = "Beans (field)",
-  vbe = "Beans (vegetables)",
-  bus = "Bushberries",
-  cab = "Cabbage",
-  car = "Carrots",
-  cit = "Citrus",
-  cot = "Cotton",
-  lin = "Linseed",
-  mai = "Maize",
-  soi = "Oil seed rape (summer)",
-  woi = "Oil seed rape (winter)",
-  oni = "Onions",
-  ape = "Peas (animals)",
-  soy = "Soybeans",
-  spr = "Spring cereals",
-  str = "Strawberries",
-  sun = "Sunflower",
-  tob = "Tobacco",
-  tom = "Tomatoes",
-  vin = "Vines")
-
-#' FOCUS PELMO crop acronyms used for naming sceneario files
-#'
-#' A named character vector with the crop acronyms used in FOCUS PELMO 5.5.3
-#' for naming the .sze files located in the FOCUS directory. The crop acronyms
-#' in the file names are sometimes capitalized, sometimes not. The scenario
-#' files used for Beans (field) and Beans (vegetable) are the same.
-#'
-#' @docType data
-#' @export
-#' @examples
-#' print(FOCUS_PELMO_crop_sze_names)
-FOCUS_PELMO_crop_sze_names = c(
-  app = "apples",
-  gra = "grass",
-  pot = "potato",
-  sug = "sbeets",
-  win = "wcerea",
-  fbe = "beans",
-  vbe = "beans", # Same sze as for fbe is used, with 'irrigation' which has no effect on the run
-  bus = "bushb",
-  cab = "cabbag",
-  car = "carrot",
-  cit = "citrus",
-  cot = "cotton",
-  lin = "linse",
-  mai = "maize",
-  soi = "rapesu",
-  woi = "rapewi",
-  oni = "onions",
-  ape = "peas",
-  soy = "soyb",
-  spr = "scerea",
-  str = "strawb",
-  sun = "sunflo",
-  tob = "tobacc",
-  tom = "tomato",
-  vin = "vines")
-
-#' Location codes in FOCUS PELMO
-#'
-#' A named character vector of one letter location codes used in FOCUS PELMO,
-#' indexed by three letter acronyms.
-
-#' @docType data
-#' @export
-#' @examples
-#' print(FOCUS_PELMO_location_codes)
-FOCUS_PELMO_location_codes = c(
-  Cha = "C", Ham = "H", Jok = "J",
-  Kre = "K", Oke = "N", Pia = "P",
-  Por = "O", Sev = "S", Thi = "T")
diff --git a/R/PELMO_runs.R b/R/PELMO_runs.R
deleted file mode 100644
index 2768205..0000000
--- a/R/PELMO_runs.R
+++ /dev/null
@@ -1,481 +0,0 @@
-#' Set up runs for FOCUS PELMO
-#'
-#' Per default, the runs are not only set up but also executed with FOCUS
-#' PELMO, the results are processed and returned. Currently, only FOCUS PELMO
-#' as installed on Linux (or other Unix systems)
-#' using the \code{install_PELMO} from the \code{PELMO.installeR} package
-#' maintained on github is supported. In such installations, FOCUS PELMO is
-#' installed into the package installation directory of \code{PELMO.installeR}
-#' and run using \code{wine}.
-#'
-#' As a side effect, an R data file (period_pfm.rda) is generated in each
-#' run directory, holding the results for all FOCUS periods, equivalent to
-#' the period.plm file generated by the FOCUS PELMO GUI.
-#'
-#' @return If evaluate is TRUE, a list of lists of matrices holding the
-#'   PEC data.
-#' @param runs A list of lists. Each inner lists has an element named 'psm'
-#'   that holds the psm string, and elements named using three letter crop acronyms,
-#'   as used in \code{\link{FOCUS_PELMO_crops}},
-#'   that hold character vectors of three letter scenario acronyms
-#'   as used in \code{\link{FOCUS_GW_scenarios_2012}}.
-#' @param psm_dir The directory where the psm files are located
-#' @param version The FOCUS PELMO version
-#' @param PELMO_base Where the FOCUS PELMO installation is located
-#' @param execute Should PELMO be executed directly?
-#' @param cores The number of cores to execute PELMO runs in parallel
-#' @param evaluate Should the results be returned?
-#' @param overwrite Should existing run directories be overwritten?
-#' @references PELMO.installeR \url{https://pkgdown.jrwb.de/PELMO.installeR}
-#'
-#' Wine \url{https://winehq.org}
-#'
-#' PELMO test results \url{http://esdac.jrc.ec.europa.eu/public_path/projects_data/focus/gw/models/pelmo/test-results/test_results_FOCUS_PELMO_5_5_3.doc}
-#' @export
-#' @examples
-#' # Reproduce the official test results for annual application of Pesticide D
-#' # to winter cereals at the day before emergence
-#' runs_1 <- list(
-#'   list(psm = 'Pesticide_D',
-#'        win = c("Cha", "Ham", "Jok", "Kre", "Oke", "Pia", "Por", "Sev", "Thi")),
-#'   list(psm = 'Pesticide_D_1_day_pre_em_every_third_year',
-#'        pot = c("Cha", "Ham", "Jok", "Kre", "Oke", "Pia", "Por", "Sev", "Thi")))
-#' time_1 <- system.time(
-#'   PECgw_1 <- PELMO_runs(runs_1, psm_dir = system.file("testdata", package = "pfm"),
-#'     cores = 15, overwrite = TRUE)
-#' )
-#' print(PECgw_1)
-#' # We get exactly the same PECgw values (on Linux, calling PELMO using Wine).
-#' print(time_1)
-#' if(!inherits(try(cpuinfo <- readLines("/proc/cpuinfo")), "try-error")) {
-#'   cat(gsub("model name\t: ", "CPU model: ", cpuinfo[grep("model name", cpuinfo)[1]]))
-#' }
-#'
-#' # Demonstrate some results with metabolites.
-#' runs_2 <- list(list(psm = 'Pesticide_D_1_May_every_other_year_mets',
-#'                     win = c("Cha", "Ham", "Kre")))
-#' PECgw_2 <- PELMO_runs(runs_2, psm_dir = system.file("testdata", package = "pfm"),
-#'   cores = 3, overwrite = TRUE)
-#' print(PECgw_2)
-PELMO_runs <- function(runs, psm_dir = ".", version = "5.5.3", PELMO_base = "auto",
-                       execute = TRUE, cores = getOption("mc.cores", 2L),
-                       evaluate = TRUE, overwrite = FALSE)
-{
-  if (PELMO_base[1] == "auto") {
-    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
-                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
-  }
-  PELMO_exe = list.files(PELMO_base, "^pelmo.*.exe", full.names = TRUE)
-  if (length(PELMO_exe) == 0) {
-    stop("Could not find PELMO executable. Did you run PELMO.installeR::install_PELMO()?")
-  }
-
-  run_list <- create_run_list(runs, psm_dir = psm_dir, check_psm_files = TRUE)
-
-  setup_run <- function(x) {
-    psm <- x$psm
-    crop <- x$crop
-    scenario <- x$scenario
-
-    location_code <- FOCUS_PELMO_location_codes[scenario]
-
-    run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
-
-    if (dir.exists(run_dir)) {
-      if (overwrite) {
-        unlink(run_dir, recursive = TRUE)
-      } else {
-        stop("Run directory for ", crop, " in ", scenario, "\n", run_dir, "\nalready exists")
-      }
-    }
-
-    # Create the run directory
-    dir.create(run_dir, recursive = TRUE)
-
-    # Copy the psm file
-    psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
-    file.copy(psm_file, run_dir)
-
-    # Copy Haude factors
-    file.copy(file.path(PELMO_base, "FOCUS", "HAUDE.DAT"), file.path(run_dir, "Haude.dat"))
-
-    # Copy scenario file
-    sze_file_upper <- paste0(location_code, "_",
-                             toupper(FOCUS_PELMO_crop_sze_names[crop]), ".sze")
-    sze_file_lower <- paste0(location_code, "_",
-                             FOCUS_PELMO_crop_sze_names[crop], ".sze")
-    if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_lower))) {
-      sze_file <- sze_file_lower
-    } else {
-      if (file.exists(file.path(PELMO_base, "FOCUS", sze_file_upper))) {
-        sze_file <- sze_file_upper
-      } else {
-        stop("Could not find szenario file for ", crop, " in ", scenario)
-      }
-    }
-    file.copy(file.path(PELMO_base, "FOCUS", sze_file), file.path(run_dir, sze_file_lower))
-
-    # Generate PELMO input file
-    input_file <- file(file.path(run_dir, "pelmo.inp"), encoding = "latin1", open = "w+")
-    on.exit(close(input_file))
-
-    add <- function(x) cat(paste0(x, "\r\n"), file = input_file, append = TRUE)
-
-    interval <- get_interval(psm_file, location_code)
-
-    n_years <- switch(as.character(interval),
-      "1" = 26,
-      "2" = 46,
-      "3" = 66)
-
-    # First line with number of climate years
-    add(paste0(formatC(n_years, width = 3, flag = "0"),
-               " 01 01 31 12 Version ", substr(version, 1, 1)))
-
-    # Second line with scenario file (why did we copy it, if this line refers to the FOCUS dir...
-    add(paste0("..\\..\\..\\", sze_file_lower))
-
-    # Third line with psm file
-    add(basename(psm_file))
-
-    # Lines with climate files
-    for (year in 1:n_years) {
-      add(paste0("..\\..\\..\\", location_code, "_", formatC(year, width = 2, flag = "0"), ".cli"))
-    }
-
-    # Output control section
-    add("00000015 00.")
-    add("    PRSN    TSER  000.0     1.0 00000001")
-    add("    TETD    TSER  000.0     1.0 00000001")
-    add("    INFL    TSER  100.0     1.0 00000001")
-    add("    RUNF    TSER  000.0     1.0 00000001")
-    add("    THET    TSER  000.0     1.0 00000001")
-    add("    THET    TSER  030.0     1.0 00000001")
-    add("    TEMP    TSER  000.0     1.0 00000001")
-    add("    TEMP    TSER  030.0     1.0 00000001")
-    add("    TPAP    TSER  000.0  1.0E05 00000001")
-    add("    TDKF    TSER  000.0  1.0E05 00000001")
-    add("    TUPF    TSER  000.0  1.0E05 00000001")
-    add("    TPST    TSER  005.0  1.0E06 00000001")
-    add("    PFLX    TSER  100.0  1.0E05 00000001")
-    add("    RFLX    TSER  000.0  1.0E05 00000001")
-    add("    LEAC    TSER  100.0  1.0E09 00000001")
-
-    # Copy pelmo executable
-    file.copy(PELMO_exe, run_dir)
-
-    # In addition to the files copied by the FOCUS PELMO GUI, we
-    # need the error file in the run directory, as we start
-    # the exe file from this directory
-    file.copy(file.path(PELMO_base, "lf90.eer"), run_dir)
-  }
-
-  lapply(run_list, setup_run)
-
-  if (execute) {
-    run_PELMO(runs, version = version, PELMO_base = PELMO_base, cores = cores)
-  }
-
-  if (evaluate) {
-    pfm_PECgw <- evaluate_PELMO(runs, version = version, PELMO_base = PELMO_base)
-    return(pfm_PECgw)
-  } else {
-    invisible(NULL)
-  }
-}
-
-#' @inheritParams PELMO_runs
-#' @importFrom parallel mclapply
-#' @rdname PELMO_runs
-#' @export
-run_PELMO <- function(runs, version = "5.5.3", PELMO_base = "auto",
-                      cores = getOption("mc.cores", 2L))
-{
-
-  if (PELMO_base[1] == "auto") {
-    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
-                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
-  }
-
-  pelmo_exe = list.files(PELMO_base, "^pelmo.*.exe")
-
-  # Create list of runs to traverse using mclapply
-  run_list <- create_run_list(runs, check_psm_files = FALSE)
-
-  # In order to run in parallel, we need to make a directory tree starting
-  # from the base of the PELMO installation, as pelmo401.exe seems to look
-  # for ..\..\..\summary.PLM and does not start if this is in use by another
-  # instance
-  run_exe <- function(x) {
-    psm <- x$psm
-    crop <- x$crop
-    scenario <- x$scenario
-    run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
-    exe_dir <- file.path(PELMO_base, paste(sample(c(0:9, letters, LETTERS),
-                                                  size = 8, replace = TRUE),
-                                           collapse = ""))
-    dir.create(exe_dir)
-
-    # Fresh PELMO base directory for every run with random name
-    run_dir_exe <- file.path(exe_dir, "FOCUS", PELMO_path(psm, crop, scenario))
-    dir.create(run_dir_exe, recursive = TRUE)
-
-    # Copy the contents of the run directory
-    run_files <- list.files(run_dir, full.names = TRUE)
-    file.copy(run_files, run_dir_exe)
-
-    # Copy FOCUS files
-    location_code <- FOCUS_PELMO_location_codes[x$scenario]
-    location_files <- list.files(file.path(PELMO_base, "FOCUS"), paste0("^", location_code),
-                                 full.names = TRUE)
-    focus_dir_exe <- file.path(exe_dir, "FOCUS")
-    file.copy(location_files, focus_dir_exe)
-    file.copy(file.path(PELMO_base, "FOCUS", "FOCUSCHECK.DAT"), focus_dir_exe)
-
-    # We need to go the directory to simplify calling pelmo with wine
-    setwd(run_dir_exe)
-    psm_file <- paste0(psm, ".psm")
-    message("Starting ", pelmo_exe, " with ", psm_file)
-    system(paste("wine", pelmo_exe, psm_file), ignore.stdout = TRUE)
-
-    # Copy the results to the original run directory
-    plm_files <- list.files(run_dir_exe, ".plm$", full.names = TRUE)
-    PLM_files <- list.files(run_dir_exe, ".PLM$", full.names = TRUE)
-    file.copy(c(plm_files, PLM_files), run_dir)
-
-    # Clean up
-    unlink(exe_dir, recursive = TRUE)
-  }
-
-  mclapply(run_list, run_exe, mc.cores = cores)
-}
-
-#' Create a path of run directories as the PELMO GUI does
-#'
-#' @export
-#' @importFrom utils data
-#' @param psm The psm identifier
-#' @param crop The PELMO crop acronym
-#' @param scenario The scenario
-PELMO_path <- function(psm, crop, scenario = NA) {
-  if (crop %in% names(FOCUS_PELMO_crops)) {
-    crop <- FOCUS_PELMO_crops[crop]
-  }
-  if (!crop %in% FOCUS_PELMO_crops) stop(crop, " is not in FOCUS_PELMO_crops.")
-  psm_dir <- paste0(psm, ".run")
-  crop_dir <- paste0(gsub(" ", "_-_", crop), ".run")
-
-  # Deal with 'irrigation'. It only affects the naming of the scenario directory,
-  # but it does not appear to affect the PELMO run. This can be seen when
-  # comparing runs set up for Beans (field) and Beans (vegetable) for the Porto
-  # scenario
-  irrigation_string <- ""
-
-  if (is.na(scenario)) {
-    return(file.path(psm_dir, crop_dir))
-  } else {
-    # 'Irrigation' is only possible in the GUI for five scenarios
-    if (scenario %in% c("Cha", "Pia", "Por", "Sev", "Thi")) {
-      crop_acronyms <- names(FOCUS_PELMO_crops)
-      names(crop_acronyms) <- FOCUS_PELMO_crops
-      crop_acronym <- crop_acronyms[[crop]]
-      # Some crops are not 'irrigated' according to the GUI
-      if (crop_acronym %in% c("win", "fbe", "woi", "ape", "spr")) {
-        irrigation_string <- "--_-_no_-_irrigation"
-      } else {
-        irrigation_string <- "--_-_irrigated"
-      }
-    }
-
-    scenario_dir <- paste0(
-      FOCUS_GW_scenarios_2012$names[scenario], "_-_(", FOCUS_PELMO_location_codes[scenario], ")",
-      irrigation_string, ".run")
-
-    return(file.path(psm_dir, crop_dir, scenario_dir))
-  }
-}
-
-#' Create a list of runs that we can traverse
-#'
-#' @inheritParams PELMO_runs
-#' @param check_psm_files Should we check if the psm file exists
-create_run_list <- function(runs, psm_dir = ".", check_psm_files = FALSE) {
-  i <- 0
-  run_list <- list()
-  for (run in runs) {
-    psm <- run$psm
-    if (check_psm_files) {
-      psm_file <- file.path(psm_dir, paste0(psm, ".psm"))
-      if (file.access(psm_file) == -1) {
-        stop(psm_file, " is not readable")
-      }
-    }
-    crops <- setdiff(names(run), "psm")
-    for (crop in crops) {
-      crop_acronyms <- names(FOCUS_PELMO_crops)
-      names(crop_acronyms) <- FOCUS_PELMO_crops
-      if (!crop %in% crop_acronyms) {
-        if (crop %in% FOCUS_PELMO_crops) {
-          crop <- crop_acronyms[crop]
-        } else {
-          stop("Invalid crop specification ", crop)
-        }
-      }
-      for (scenario in run[[crop]]) {
-        i <- i + 1
-        run_list[[i]] <- list(psm = psm, crop = crop, scenario = scenario)
-      }
-    }
-  }
-  return(run_list)
-}
-
-#' @rdname PELMO_runs
-#' @inheritParams PELMO_runs
-#' @importFrom parallel mclapply
-#' @export
-evaluate_PELMO <- function(runs, version = "5.5.3", PELMO_base = "auto")
-{
-
-  if (PELMO_base[1] == "auto") {
-    PELMO_base = file.path(system.file(package = "PELMO.installeR"),
-                                       paste0("FOCUSPELMO.", gsub("\\.", "", version)))
-  }
-
-  pfm_PECgw <- list()
-  for (run in runs) {
-    psm <- run$psm
-    pfm_PECgw[[psm]] <- list()
-    crops <- setdiff(names(run), "psm")
-
-    # Get acronyms of simulated compounds
-    example_run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crops[1], run[[crops[1]]][1]))
-    example_echo_file <- readLines(file.path(example_run_dir, "ECHO.PLM"), encoding = "latin1")
-    parm_lines <- grep("\\*\\*\\* PARAMETERS OF", example_echo_file, value = TRUE)
-    acronyms <- gsub(".*\\((.*)\\).*", "\\1", parm_lines)
-    met_codes <- gsub(".*METABOLITE (..).*", "\\1", parm_lines)
-    met_codes[1] <- NA
-    names(met_codes) <- acronyms
-
-    # Loop over runs to get results
-    for (crop in crops) {
-      scenarios <- run[[crop]]
-
-      pfm_PECgw[[psm]][[crop]] <- matrix(nrow = length(scenarios), ncol = length(acronyms),
-                                         dimnames = list(scenarios, acronyms))
-
-      for (scenario in scenarios) {
-        run_dir <- file.path(PELMO_base, "FOCUS", PELMO_path(psm, crop, scenario))
-
-        psm_file <- file.path(run_dir, paste0(psm, ".psm"))
-        location_code <- FOCUS_PELMO_location_codes[scenario]
-        interval <- get_interval(psm_file, location_code)
-
-        # Get percolate for each period
-        wasser <- readLines(file.path(run_dir, "WASSER.PLM"))
-        percolate_lines <- grep("RECHARGE BELOW ROOT ZONE", wasser, value = TRUE)
-        percolate <- as.numeric(substr(percolate_lines, 40, 46))
-        percolate_period <- sum_periods(percolate, interval)
-
-        # Set up results that should match period.plm generated by the PELMO GUI
-        results_pfm <- list()
-        for (acronym in acronyms) {
-          results_pfm[[acronym]] <- list()
-          periods <- data.frame(
-            period = 1:20,
-            flux = NA,
-            percolate = 10 * percolate_period,
-            conc = NA)
-
-          if (is.na(met_codes[acronym])) {
-            chem_file <- file.path(run_dir, "CHEM.PLM")
-          } else {
-            chem_file <- file.path(run_dir, paste0("CHEM_", met_codes[acronym], ".PLM"))
-          }
-
-          annual_flux <- get_flux(chem_file)
-          periods$flux <- sum_periods(annual_flux, interval) * 1000
-
-          periods$conc <- 100 * periods$flux / periods$percolate
-
-          results_pfm[[acronym]]$periods <- periods
-          PECgw <- focus_80th(periods$conc)
-
-          results_pfm[[acronym]]$focus <- PECgw
-
-          pfm_PECgw[[psm]][[crop]][scenario, acronym] <- round(PECgw, 3)
-        }
-        save(results_pfm, file = file.path(run_dir, "period_pfm.rda"))
-      }
-    }
-  }
-  return(pfm_PECgw)
-}
-
-#' Get the application interval in years from a psm file
-#'
-#' @param psm_file The path to the .psm file
-#' @param location_code The location code
-get_interval <- function(psm_file, location_code) {
-  # How many years do we calculate (26, 46 or 66)?
-  psm <- readLines(psm_file, encoding = "latin1")
-  number_of_apps_lines <- grep("number of application location", psm)
-  absolute_apps_line <- grep(location_code,
-                             psm[number_of_apps_lines])
-  interval <- 1 # application every year
-  if (length(absolute_apps_line) == 1) {
-    apps_root <- number_of_apps_lines[absolute_apps_line]
-  } else {
-    apps_root <- number_of_apps_lines[1]
-  }
-
-  number_of_apps <- as.integer(substr(psm[apps_root], 1, 3))
-  last_app_line <- psm[apps_root + number_of_apps]
-  last_app_year <- as.integer(gsub("^.{2,3}  ..  (..)   .*", "\\1",
-                                   last_app_line))
-  if (last_app_year > 26) interval <- 2
-  if (last_app_year > 46) interval <- 3
-  return(interval)
-}
-
-#' Sum up values according to FOCUS periods
-#' 
-#' @param annual The annual flux as obtained by \code{get_flux}
-#' @param interval The interval in years
-sum_periods <- function(annual, interval) {
-  n_years <- switch(as.character(interval),
-    "1" = 26,
-    "2" = 46,
-    "3" = 66)
-  period_start_years <- seq(from = 7, to = n_years, by = interval)
-  sapply(1:20, function(x) {
-    years_i <- period_start_years[x] + (0 : (interval - 1))
-    sum(annual[years_i])
-  })
-}
-
-#' Get the flux of a chemical out of the FOCUS layer from a CHEM*.PLM file
-#' 
-#' @param chem_file The full path to a CHEM*.PLM file
-get_flux <- function(chem_file) {
-  chem <- readLines(chem_file)
-  lowest_focus_comp_lines <- grep("^    .     21      ", chem, value = TRUE)
-  lowest_focus_comp <- read.table(text = lowest_focus_comp_lines)
-  return(lowest_focus_comp$V9)
-}
-
-#' Calculate the 80th percentile according to FOCUS guidance
-#'
-#' This is nowadays defined as the mean of the 16th and the 17th
-#' highest value. Previously, the 17th highest values was used (FOCUS 2014, p.
-#' 18). NaN values need to be set to zero in order to reproduce the
-#' values obtained by PELMO.
-#' @param c_period A numeric vector of values to calculate the percentile from
-#' @param old Should the old calculation method be used (the 17th highest value)?
-#' @export
-focus_80th <- function(c_period, old = FALSE) {
-  c_period <- ifelse(is.na(c_period), 0, c_period)
-   c_period_sorted <- sort(c_period)
-  if (old) return(c_period_sorted[17])
-  else return(mean(c_period_sorted[c(16, 17)]))
-}