R/util.units.R

In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

# CHNOSZ/util.units.R
# Set units and convert values between units

P.units <- function(units = NULL) {
  ## Change units of pressure or list the current one
  # Show the current units, if none is specified
  if(is.null(units)) return(get("thermo", CHNOSZ)$opt$P.units)
  # Argument handling
  units <- tolower(units)
  if(!units %in% c("bar", "mpa")) stop("units of pressure must be either bar or MPa")
  # Set the units and return them
  if(units == "bar") with(CHNOSZ, thermo$opt$P.units <- "bar")
  if(units == "mpa") with(CHNOSZ, thermo$opt$P.units <- "MPa")
  message("changed pressure units to ", get("thermo", CHNOSZ)$opt$P.units)
}

T.units <- function(units = NULL) {
  ## Change units of temperature or list the current one
  # Show the current units, if none is specified
  if(is.null(units)) return(get("thermo", CHNOSZ)$opt$T.units)
  # Argument handling
  units <- tolower(units)
  if(!units %in% c("c", "k")) stop("units of temperature must be either C or K")
  # Set the units and return them
  if(units == "c") with(CHNOSZ, thermo$opt$T.units <- "C")
  if(units == "k") with(CHNOSZ, thermo$opt$T.units <- "K")
  message("changed temperature units to ", get("thermo", CHNOSZ)$opt$T.units)
}

E.units <- function(units = NULL) {
  ## Change units of energy or list the current one
  # Show the current units, if none is specified
  if(is.null(units)) return(get("thermo", CHNOSZ)$opt$E.units)
  # Argument handling
  units <- tolower(units)
  if(!units %in% c("cal", "j")) stop("units of energy must be either cal or J")
  # Set the units and return them
  if(units == "cal") with(CHNOSZ, thermo$opt$E.units <- "cal")
  if(units == "j") with(CHNOSZ, thermo$opt$E.units <- "J")
  message("changed energy units to ", get("thermo", CHNOSZ)$opt$E.units)
}

convert <- function(value, units, T = 298.15,
  P = 1, pH = 7, logaH2O = 0) {
  # Converts value(s) to the specified units

  # Process a list value if it's the output from solubility 20190525
  if(is.list(value) & !is.data.frame(value)) {
    if(!isTRUE(value$fun %in% c("solubility", "solubilities"))) stop("'value' is a list but is not the output from solubility()")
    if(!is.character(units)) stop("please specify a character argument for the destination units (e.g. ppm or logppm)")
    # Determine the element from 'balance' or 'in.terms.of', if it's available
    element <- value$in.terms.of
    if(is.null(element)) element <- value$balance
    grams.per.mole <- mass(element)
    message(paste("solubility: converting to", units, "by weight using the mass of", element))
    ppfun <- function(loga, units, grams.per.mole) {
      # Exponentiate loga to get molality
      moles <- 10^loga
      # Convert moles to mass (g)
      grams <- moles * grams.per.mole
      # Convert grams to ppb, ppm, or ppt
      ppx <- NULL
      # 1 ppt = 1 g / kg H2O
      # 1 ppm = 1 mg / kg H2O
      if(grepl("ppt", units)) ppx <- grams * 1e0
      if(grepl("ppm", units)) ppx <- grams * 1e3
      if(grepl("ppb", units)) ppx <- grams * 1e6
      if(is.null(ppx)) stop(paste("units", units, "not available for conversion"))
      # Use the logarithm if specified
      if(grepl("log", units)) ppx <- log10(ppx)
      ppx
    }
    # Do the conversion for the conserved basis species, then each species
    value$loga.balance <- ppfun(value$loga.balance, units, grams.per.mole)
    value$loga.equil <- lapply(value$loga.equil, ppfun, units = units, grams.per.mole = grams.per.mole)
    # Identify the units in the function text
    value$fun <- paste(value$fun, units, sep = ".")
    # Return the updated object
    return(value)
  }

  ### Argument handling for non-list value
  if(is.null(value)) return(NULL)
  if(!is.character(units)) stop(paste('convert: please specify',
    'a character argument for the destination units.\n',
    'possibilities include (G or logK) (C or K) (J or cal) (cm3bar or calories) (Eh or pe)\n',
    'or their lowercase equivalents.\n'), call. = FALSE)
  Units <- units # For the possible message to user
  units <- tolower(units)

  # Tests and calculations for the specified units
  if(units %in% c('c', 'k')) {
    CK <- 273.15
    if(units == 'k') value <- value + CK
    if(units == 'c') value <- value - CK 
  }
  else if(units[1] %in% c('j', 'cal')) {
    Jcal <- 4.184
    if(units == 'j') value <- value * Jcal
    if(units == 'cal') value <- value / Jcal
  }
  else if(units %in% c('g', 'logk')) {
    # Gas constant
    #R <- 1.9872  # cal K^-1 mol^-1  Value used in SUPCRT92
    #R <- 8.314445  # = 1.9872 * 4.184 J K^-1 mol^-1  20220325
    R <- 8.314463  # https://physics.nist.gov/cgi-bin/cuu/Value?r 20230630
    if(units == 'logk') value <- value / (-log(10) * R * T)
    if(units == 'g') value <- value * (-log(10) * R * T)
  }
  else if(units %in% c('cm3bar', 'joules')) {
    if(units == 'cm3bar') value <- value * 10
    if(units == 'joules') value <- value / 10
  }
  else if(units %in% c('eh', 'pe')) {
    R <- 0.00831470
    F <- 96.4935
    if(units == 'pe') value <- value * F / ( log(10) * R * T )
    if(units == 'eh') value <- value * ( log(10) * R * T ) / F
  }
  else if(units %in% c('bar', 'mpa')) {
    barmpa <- 10
    if(units == 'mpa') value <- value / barmpa
    if(units == 'bar') value <- value * barmpa
  }
  else if(units %in% c('e0', 'logfo2')) {
    # Convert between Eh and logfO2
    supcrt.out <- suppressMessages(subcrt(c("H2O", "oxygen", "H+", "e-"), c(-1, 0.5, 2, 2), T = T, P = P, convert = FALSE))
    if(units == 'logfo2') value <- 2*(supcrt.out$out$logK + logaH2O + 2*pH + 2*(convert(value, 'pe', T = T)))
    if(units == 'e0') value <- convert(( -supcrt.out$out$logK - 2*pH + value/2 - logaH2O )/2, 'Eh', T = T)
  }
  else cat(paste('convert: no conversion to ', Units, ' found.\n', sep = ''))
  return(value)
}

### Unexported functions ###

outvert <- function(value, units) {
  # Converts the given value from the given units to those specified in thermo()$opt
  units <- tolower(units)
  opt <- get("thermo", CHNOSZ)$opt
  if(units %in% c("c", "k")) {
    if(units == "c" & opt$T.units == "K") return(convert(value, "k"))
    if(units == "k" & opt$T.units == "C") return(convert(value, "c"))
  }
  if(units %in% c("j", "cal")) {
    if(units == "j" & opt$E.units == "cal") return(convert(value, "cal"))
    if(units == "cal" & opt$E.units == "J") return(convert(value, "j"))
  }
  if(units %in% c("bar", "mpa")) {
    if(units == "mpa" & opt$P.units == "bar") return(convert(value, "bar"))
    if(units == "bar" & opt$P.units == "MPa") return(convert(value, "mpa"))
  }
  return(value)
}

envert <- function(value, units) {
  # Convert values to the specified units
  # from those given in thermo()$opt
  if(!is.numeric(value[1])) return(value)
  units <- tolower(units)
  opt <- get("thermo", CHNOSZ)$opt
  if(units %in% c('c', 'k', 't.units')) {
    if(units == 'c' & opt$T.units == 'K') return(convert(value, 'c'))
    if(units == 'k' & opt$T.units == 'C') return(convert(value, 'k'))
  }
  if(units %in% c('j', 'cal', 'e.units')) {
    if(units == 'j' & opt$T.units == 'Cal') return(convert(value, 'j'))
    if(units == 'cal' & opt$T.units == 'J') return(convert(value, 'cal'))
  }
  if(units %in% c('bar', 'mpa', 'p.units')) {
    if(units == 'mpa' & opt$P.units == 'bar') return(convert(value, 'mpa'))
    if(units == 'bar' & opt$P.units == 'MPa') return(convert(value, 'bar'))
  }
  return(value)
}