R/calc_OSLLxTxDecomposed.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Calculate Lx/Tx ratio for decomposed CW-OSL signal components
#'
#' @description Calculate `Lx/Tx` ratios from a given set of decomposed
#' CW-OSL curves decomposed by `[OSLdecomposition::RLum.OSL_decomposition]`
#'
#' @param OSL.component [integer] or [character] (*optional*):
#' a single index or a name describing which OSL signal component shall be evaluated.
#' This argument can either be the name of the OSL component assigned by
#' `[OSLdecomposition::RLum.OSL_global_fitting]` or the index of component.
#' Then `'1'` selects the fastest decaying component, `'2'` the
#' second fastest and so on. If not defined, the fastest decaying component is selected.
#'
#' @param Lx.data [data.frame] (**required**): Component table created by
#' `[OSLdecomposition::RLum.OSL_decomposition]` and per default located
#' at `object@records[[...]]@info$COMPONENTS`.The value of `$n[OSL.component]`
#' is set as `LnLx`. The value of `$n.error[OSL.component]` is set as `LnLx.error`
#'
#' @param Tx.data [data.frame] (*optional*): Component table created by
#' `[OSLdecomposition::RLum.OSL_decomposition]` and per default located at
#' `object@records[[...]]@info$COMPONENTS`. The value of `$n[OSL.component]`
#' is set as `TnTx`. The value of `$n.error[OSL.component]` is set as `TnTx.error`
#'
#' @param sig0 [numeric] (*with default*): allows adding an extra error component
#' to the final `Lx/Tx` error value (e.g., instrumental error).
#'
#' @param digits [integer] (*with default*): round numbers to the specified digits.
#' If digits is set to `NULL` nothing is rounded.
#'
#' @return Returns an S4 object of type [RLum.Results-class].
#'
#' Slot `data` contains a [list] with the following structure:
#'
#' **@data**
#' ```
#' $LxTx.table (data.frame)
#' .. $ LnLx
#' .. $ TnTx
#' .. $ Net_LnLx
#' .. $ Net_LnLx.Error
#' .. $ Net_TnTx
#' .. $ Net_TnTx.Error
#' .. $ LxTx
#' .. $ LxTx.relError
#' .. $ LxTx.Error
#' ```
#'
#' @section Function version: 0.1.0
#'
#' @author Dirk Mittelstrass
#'
#' @seealso [RLum.Data.Curve-class], [plot_GrowthCurve], [analyse_SAR.CWOSL]
#'
#' @references Mittelstrass D., Schmidt C., Beyer J., Straessner A., 2019.
#' Automated identification and separation of quartz CW-OSL signal components with R.
#' talk presented at DLED 2019, Bingen, Germany
#' [http://luminescence.de/OSLdecomp_talk.pdf]()\cr
#'
#' @keywords datagen
#' @md
#' @export
calc_OSLLxTxDecomposed <- function(
  Lx.data,
  Tx.data = NULL,
  OSL.component = 1L,
  sig0 = 0,
  digits = NULL
){

  # ToDo:
  # - Integrity checks for the component table
  # - Handle background-signal-component if present
  # - add Tx.data integrity checks
  # - add previous-residual-subtraction functionality
  # - add list with decomposition algorithm parameters to return object
  # - add example in documentation

  ##--------------------------------------------------------------------------##
  ## (1) - integrity checks
  if (!(is.data.frame(Lx.data) && (nrow(Lx.data) >= 1)))
    stop("[calc_OSLLxTxDecomposed()] No valid component data.frame for Lx value", call. = FALSE)

  if (!(is.null(Tx.data)) && !(is.data.frame(Tx.data) && (nrow(Tx.data) >= 1)))
    stop("[calc_OSLLxTxDecomposed()] No valid component data.frame for Tx value", call. = FALSE)

  # define the component
  component_index <- NA

  #select only the first element; we do this silently because it is clearly
  #written in the documentation
  OSL.component <- as.integer(OSL.component[1])

  if (!(is.numeric(OSL.component) || is.character(OSL.component)))
    stop("[calc_OSLLxTxDecomposed()] Type error! No valid data type for OSL.component", call. = FALSE)

  # get component index from component name
  if (is.character(OSL.component)) {
    if (tolower(OSL.component) %in% tolower(Lx.data$name)) {
      component_index <- which(tolower(OSL.component) == tolower(Lx.data$name))

    } else {
      stop(paste0("[calc_OSLLxTxDecomposed()] Invalid OSL component name! Valid names are: ",
                  paste(Lx.data$name, collapse = ", ")), call. = FALSE)
    }
  }

  # if a numeric is given, check if it matches with any component index
  if (is.numeric(OSL.component)) {
    if (OSL.component %in% 1:nrow(Lx.data)) {
      component_index <- OSL.component

      # insert background-signal-component check here

    } else {
      stop(paste0("[calc_OSLLxTxDecomposed()] Invalid OSL component index!
                  Component table has ", nrow(Lx.data), " rows."))

    }

  }

  ##--------------------------------------------------------------------------##
  ## (2) - extract Lx and Tx values

  LnLx <- Lx.data$n[component_index]
  LnLx.Error <- Lx.data$n.error[component_index]

  TnTx <- 1
  TnTx.Error <- 0
  if (!is.null(Tx.data)) {

    TnTx <- Tx.data$n[component_index]
    TnTx.Error <- Tx.data$n.error[component_index]
  }

  ##combine results
  LnLxTnTx <- cbind(
    LnLx,
    LnLx.Error,
    TnTx,
    TnTx.Error
  )

  # THE FOLLOWING CODE IS MOSTLY IDENTICAL WITH (4) IN calc_OSLLxTxRatio()

  ##--------------------------------------------------------------------------##
  ##(4) Calculate LxTx error according Galbraith (2014)

  #transform results in a data.frame
  LnLxTnTx <- as.data.frame(LnLxTnTx)

  #add col names
  colnames(LnLxTnTx)<-c("Net_LnLx", "Net_LnLx.Error",
                        "Net_TnTx", "Net_TnTx.Error")

  ##calculate Ln/Tx
  LxTx <- LnLxTnTx$Net_LnLx/LnLxTnTx$Net_TnTx

  ##set NaN
  if(is.nan(LxTx)) LxTx <- 0

  ##calculate Ln/Tx error
  LxTx.relError <- sqrt((LnLx.Error / LnLx)^2 + (TnTx.Error / TnTx)^2)
  LxTx.Error <- abs(LxTx * LxTx.relError)

  ##set NaN
  if(is.nan(LxTx.Error)) LxTx.Error <- 0

  ##add an extra component of error
  LxTx.Error <- sqrt(LxTx.Error^2 + (sig0 * LxTx)^2)

  ##return combined values
  temp <- cbind(LnLxTnTx, LxTx, LxTx.Error)


  ##apply digits if wanted
  if(!is.null(digits)){
    temp[1,] <- round(temp[1,], digits = digits)

  }

  # ToDo: Add decomposition algorithm parameters here
  # calc.parameters <- list(...)

  ##set results object
  return(set_RLum(
      class = "RLum.Results",
      data = list(
        LxTx.table = temp),
      #  calc.parameters = calc.parameters),
      info = list(call = sys.call())
  ))

}