R function to calculate equivalent black carbon (EBC) components with the aethalometer model.

calculate_ebc_components.R

#' Function to calculate equivalent black carbon (EBC) components with the 
#' aethalometer model.
#' 
#' @param df Data frame/tibble containing aethalometer absorption observations. 
#' See examples for details of format required. 
#' 
#' @param uv Name of variable/column in \code{df} to be used for the UV 
#' absorption.
#' 
#' @param ir Name of variable/column in \code{df} to be used for the IR 
#' absorption.
#' 
#' @param slope A model slope from an absorption ~ elemental carbon linear 
#' regression model. This is also called the mass absorption cross section/mass 
#' absorption coefficient (MAC) value. This value can either be a single value,
#' or a possibly dynamic value in \code{df}. If the variable exists in \code{df},
#' use the variable name. 
#' 
#' @param uv_nm UV wavelength, this should match \code{uv}.
#' 
#' @param ir_nm IR wavelength, this should match \code{ir}.
#' 
#' @param alpha_tr Angstrom exponent (usually denoted by alpha) used to 
#' determine the traffic contribution of EBC. 
#' 
#' @param alpha_wb Angstrom exponent (usually denoted by alpha) used to 
#' determine the woodburnng contribution of EBC. 
#' 
#' @author Stuart K. Grange
#' 
#' @return Tibble with three variables, EBC WB (woodburning), EBC TR (traffic), 
#' and EBC total mass concentrations. 
#' 
#' @examples 
#' 
#' # Example data frame
#' data_example <- structure(
#'   list(
#'     date = structure(
#' c(1546300800L, 1546304400L, 1546308000L, 1546311600L, 1546315200L), 
#' class = c("POSIXct", "POSIXt"), 
#' tzone = "UTC"
#' ), 
#' absorption_470_nm = c(51.74, 46.5, NA, 23.7675, 8.73),
#' absorption_950_nm = c(20.39, 19.59, NA, 7.135, 3.22)),
#' class = c("tbl_df", "tbl", "data.frame" ), 
#'   row.names = c(NA, -5L)
#' )
#' 
#' # Use function
#' calculate_ebc_components(
#'   data_example, 
#'   uv = "absorption_470_nm", 
#'   ir = "absorption_950_nm",
#'   slope = 10,
#'   uv_nm = 470, 
#'   ir_nm = 950, 
#'   alpha_tr = 0.9, 
#'   alpha_wb = 1.68
#' )
#' 
calculate_ebc_components <- function(df, uv = "absorption_470_nm", 
                                     ir = "absorption_950_nm",
                                     slope, uv_nm = 470, ir_nm = 950, 
                                     alpha_tr = 0.9, alpha_wb = 1.68) {
  
  # Get values as vector
  uv_values <- df %>% 
    select(!!uv) %>% 
    pull()
  
  ir_values <- df %>% 
    select(!!ir) %>% 
    pull()
  
  # Get vector of slopes from data frame, if a single value was not supplied
  if (!class(slope) == "numeric" && length(slope) == 1) {
    slope <- df %>% 
      select(!!slope) %>% 
      pull()
  }
  
  # Calculate relations, used a few times later
  relation_tr <- (uv_nm / ir_nm) ^ -alpha_tr
  relation_wb <- (uv_nm / ir_nm) ^ -alpha_wb
  
  # Calculate absorptions
  adsorption_ir_tr <- (uv_values - ir_values * relation_wb) / 
    (relation_tr - relation_wb)
  
  adsorption_ir_wb <- (uv_values - ir_values * relation_tr) / 
    (relation_wb - relation_tr)

  # Calculate ebc components
  ebc_tr <- adsorption_ir_tr / slope
  ebc_wb <- adsorption_ir_wb / slope
  ebc_total <- ebc_wb + ebc_tr 
  
  # Build tibble for return
  df <- tibble(ebc_tr, ebc_wb, ebc_total)
  
  return(df)
  
}