R/calculate_rg.R

In bentrueman/fffprocessr: Process FFF-UV-MALS-ICP-MS Data

#' Calculate the radius of gyration from FFF-MALS data
#'
#' @param data A tibble generated by `load_mals()`.
#' @param window The time interval on which to return calculated radii.
#' @param method Use method "zimm" or "watt". Method "watt" is experimental.
#' @param calibration A data frame with three columns: number (angle number), theta (angle), and
#' coefficient (MALS instrument calibration constant).
#' @param order Order of polynomial fit for method = "watt"
#'
#' @return A tibble generated by `load_mals()` with additional columns 'timeslice', 'theta',
#' 'rayleigh_ratio', and either 'rg_zimm' or 'rg_watt'.
#' @export
#' @importFrom dplyr %>%
#' @importFrom rlang .data
#' @examples
#' path <- system.file("extdata/mals", package = "fffprocessr")
#' mals_data <- load_mals(path = path)
#' calculate_rg(mals_data)
calculate_rg <- function(
  data,
  window = .1,
  method = "zimm",
  order = 2,
  calibration = mals_calib
) {
  # mals calibration coefficients:
  mals_calib <- calibration

  data %>%
    dplyr::mutate(
      theta = stringr::str_extract(.data$param, "\\d+") %>% as.numeric(),
      theta_rad = pi * .data$theta / 180, # convert to radians
      x = sin(.data$theta_rad / 2) ^ 2
    ) %>%
    # MALS calibration coefficients
    dplyr::left_join(mals_calib, by = "theta") %>%
    dplyr::mutate(
      rayleigh_ratio = .data$conc * .data$coefficient,
      y = 1 / .data$rayleigh_ratio
    ) %>%
    dplyr::group_by(date, sample, timeslice = plyr::round_any(.data$time, window)) %>%
    tidyr::nest() %>%
    dplyr::ungroup() %>%
    dplyr::mutate(
      model = purrr::map(data, ~ if(method == "zimm") {stats::lm(y ~ x, data = .x)} else {
        stats::lm(rayleigh_ratio ~ poly(x, order, raw = TRUE), data = .x)
      }),
      coef = purrr::map(.data$model, ~ stats::coef(.x)),
      # include the index of refraction of water (1.33)?
      # Many texts leave it out, but it does result in a match to the PostNova program output
      rg_zimm = purrr::map(
        .data$coef,
        ~ if(method == "zimm") {
          suppressWarnings(sqrt(3 * 532 ^ 2 * (.x[2] / .x[1]) / (16 * pi ^ 2 * 1.33 ^ 2)))
        } else NA_real_
      ),
      # from https://doi.org/10.1007/s11051-018-4397-x
      rg_watt = purrr::map(
        .data$coef,
        ~ if(method == "watt") {
          suppressWarnings(sqrt(3 * 532 ^ 2 * (-.x[2]/.x[1])/ (16 * pi ^ 2)))
        } else NA_real_
      ),
    ) %>%
    tidyr::unnest(c(tidyselect::starts_with("rg"), data)) %>%
    dplyr::select_if(~ !is.list(.x)) %>%
    rg_select(method)
}

rg_select <- function(x, method) {
  param <- NULL
  timeslice <- NULL
  time <- NULL
  conc <- NULL
  theta <- NULL
  rayleigh_ratio <- NULL
  rg_zimm <- NULL
  rg_watt <- NULL
  if(method == "zimm") {
    dplyr::select(x,
                  file, date, sample, param, timeslice, time, conc, theta,
                  rayleigh_ratio, rg_zimm
    )
  } else if(method == "watt") {
    dplyr::select(x,
                  file, date, sample, param, timeslice, time, conc, theta,
                  rayleigh_ratio, rg_watt
    )
  } else "Choose either method 'zimm' or 'watt'"
}