R/EOSregress.R

In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

# CHNOSZ/EOSregress.R  
# Model volumes and heat capacities of aqueous species
# 20091105 first version
# 20110429 revise and merge with CHNOSZ package

Cp_s_var <- function(T = 298.15, P = 1, omega.PrTr = 0, Z = 0) {
  # Solvation contribution to heat capacity in the HKF EOS, divided by omega(Pr,Tr) (Joules)
  Cp_s <- hkf("Cp", parameters = data.frame(omega = omega.PrTr, Z = Z), T = T, P = P, contrib = "s")$aq
  return(Cp_s[[1]][, 1] / omega.PrTr)
}

V_s_var <- function(T = 298.15, P = 1, omega.PrTr = 0, Z = 0) {
  # Solvation contribution to volume in the HKF EOS, divided by omega(Pr,Tr) (cm3.bar)
  # [the negative sign on this term as written in the HKF EOS is accounted for by hkf()]
  V_s <- hkf("V", parameters = data.frame(omega = omega.PrTr, Z = Z), T = T, P = P, contrib = "s")$aq
  return(V_s[[1]][, 1]/convert(omega.PrTr, "cm3bar"))
}

EOSvar <- function(var, T, P, ...) {
  # Get the variables of a term in a regression equation
  # T (K), P (bar)
  Theta <- 228 # K
  Psi <- 2600  # bar
  out <- switch(EXPR = var,
    "(Intercept)" = rep(1, length(T)),
    "T" = T,
    "P" = P,
    "TTheta" = T - Theta,                 # T-Theta
    "invTTheta" = (T - Theta)^-1,         # 1/(T-Theta)
    "TTheta2" = (T - Theta)^2,            # (T-Theta)^2
    "invTTheta2" = (T - Theta)^-2,        # 1/(T-Theta)^2
    "invPPsi" = (P + Psi)^-1,             # 1/(P+Psi)
    "invPPsiTTheta" = (P + Psi)^-1 * (T - Theta)^-1,  # 1/[(P+Psi)(T-Theta)]
    "TXBorn" = T*water("XBorn", T = T, P = P)[, 1],
    "drho.dT" = -water("rho", T = T, P = P)[, 1]*water("E", T = T, P = P)[, 1],
    "V.kT" = water("V", T = T, P = P)[, 1]*water("kT", T = T, P = P)[, 1],
    # Fallback: get a variable that is a property of water, or
    # is any other function by name (possibly a user-defined function)
    (  if(var %in% water.SUPCRT92()) water(var, T, P)[, 1]
       else if(exists(var)) {
         if(is.function(get(var))) {
           if(all(c("T", "P") %in% names(formals(get(var))))) get(var)(T = T, P = P, ...)
           else stop(paste("the arguments of ", var, "() do not contain T and P", sep = ""))
         }
         else stop(paste("an object named", var, "is not a function"))
       }
       else stop(paste("can't find a variable named", var))
    )
  )
  # 20151126 Apply the negative sign in the HKF EOS for V to the variable
  # (not to omega as previously assumed)
  if(var == "QBorn") out <- -out
  return(out)
}

EOSlab <- function(var, coeff = "") {
  # Make pretty labels for the variables
  lab <- switch(EXPR = var,
    # These are regression variables listed in EOSregress.Rd
    "(Intercept)" = substitute(YYY*" ", list(YYY = coeff)),
    "T" = substitute(YYY%*%italic(T), list(YYY = coeff)),
    "P" = substitute(YYY%*%italic(P), list(YYY = coeff)),
    "TTheta" = substitute(YYY%*%(italic(T)-Theta), list(YYY = coeff)),
    "invTTheta" = substitute(YYY/(italic(T)-Theta), list(YYY = coeff)),
    "TTheta2" = substitute(YYY%*%(italic(T)-Theta)^2, list(YYY = coeff)),
    "invTTheta2" = substitute(YYY/(italic(T)-Theta)^2, list(YYY = coeff)),
    "invPPsi" = substitute(YYY/(italic(P)+Psi),list(YYY = coeff)),
    "invPPsiTTheta" = substitute(YYY/((italic(P)+Psi)(italic(T)-Theta)), list(YYY = coeff)),
    "TXBorn" = substitute(YYY%*%italic(TX), list(YYY = coeff)),
    "drho.dT" = substitute(YYY%*%(d~rho/dT), list(YYY = coeff)),
    "V.kT" = substitute(YYY%*%V~kappa[italic(T)], list(YYY = coeff)),
    # These are non-single-letter properties of water as listed in water.Rd
    "kT" = substitute(YYY%*%kappa[italic(T)], list(YYY = coeff)),
    "alpha" = substitute(YYY%*%alpha, list(YYY = coeff)),
    "beta" = substitute(YYY%*%beta, list(YYY = coeff)),
    "epsilon" = substitute(YYY%*%epsilon, list(YYY = coeff)),
    "rho" = substitute(YYY%*%rho, list(YYY = coeff)),
    "NBorn" = substitute(YYY%*%italic(N), list(YYY = coeff)),
    "QBorn" = substitute(YYY%*%italic(Q), list(YYY = coeff)),
    "XBorn" = substitute(YYY%*%italic(X), list(YYY = coeff)),
    "YBorn" = substitute(YYY%*%italic(Y), list(YYY = coeff)),
    "ZBorn" = substitute(YYY%*%italic(Z), list(YYY = coeff)),
    (
      # If var is a function, does have an attribute named "label"?
      if(exists(var)) {
        if(is.function(get(var))) {
          if(!is.null(attr(get(var), "label"))) {
            return(substitute(YYY*XXX, list(YYY = coeff, XXX = attr(get(var), "label"))))
            # Fallback, use the name of the variable
            # (e.g. for a property of water such as A, G, S, U, H, or name of a user-defined function)
          } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
        } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
      } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
    )
  )
  return(lab)
}

EOSregress <- function(exptdata, var = "", T.max = 9999, ...) {
  # Regress exptdata using terms listed in fun 
  # Which values to use
  iT <- which(exptdata$T <= T.max)
  exptdata <- exptdata[iT, ]
  # Temperature and pressure
  T <- exptdata$T
  P <- exptdata$P
  # The third column is the property of interest: Cp or V
  X <- exptdata[, 3]
  # Now build a regression formula 
  if(length(var) == 0) stop("var is missing")
  fmla <- as.formula(paste("X ~ ", paste(var, collapse = "+")))
  # Retrieve the values of the variables
  for(i in seq_along(var)) assign(var[i], EOSvar(var[i], T = T, P = P, ...))
  # Now regress away!
  EOSlm <- lm(fmla)
  return(EOSlm)
}

EOScalc <- function(coefficients, T, P, ...) {
  # Calculate values of volume or heat capacity from regression fit
  X <- 0
  for(i in 1:length(coefficients)) {
    coeff.i <- coefficients[[i]]
    fun.i <- EOSvar(names(coefficients)[i], T, P, ...)
    X <- X + coeff.i * fun.i
  }
  return(X)
}

EOSplot <- function(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
  fun.legend = "topleft", coefficients = NULL, add = FALSE,
  lty = par("lty"), col = par("col"), ...) {
  # Plot experimental and modelled volumes and heat capacities
  # First figure out the property (Cp or V) from the exptdata
  prop <- colnames(exptdata)[3]
  # If var is NULL use HKF equations
  if(is.null(var)) {
    if(prop == "Cp") var <- c("invTTheta2","TXBorn")
    if(prop == "V") var <- c("invTTheta","QBorn")
  }
  # Perform the regression, only using temperatures up to T.max
  if(is.null(coefficients)) {
    EOSlm <- EOSregress(exptdata, var, T.max, ...)
    coefficients <- EOSlm$coefficients
  }
  # Only plot points below a certain temperature
  iexpt <- 1:nrow(exptdata)
  if(!is.null(T.plot)) iexpt <- which(exptdata$T < T.plot)
  # For a nicer plot, extend the ranges, but don't go below -20 degrees C
  ylim <- extendrange(exptdata[iexpt, prop], f = 0.1)
  xlim <- extendrange(exptdata$T[iexpt], f = 0.1)
  xlim[xlim < 253.15] <- 253.15
  # Start plot
  if(!add) {
    thermo.plot.new(xlim = xlim, ylim = ylim, xlab = axis.label("T", units = "K"),
      ylab = axis.label(paste(prop, "0", sep = "")), yline = 2, mar = NULL)
    # Different plot symbols to represent size of residuals
    pch.open <- 1
    pch.filled <- 16
    # Find the calculated values at these conditions
    calc.X <- EOScalc(coefficients, exptdata$T, exptdata$P, ...)
    expt.X <- exptdata[, prop]
    # Are we within 10% of the values?
    in10 <- which(abs((calc.X-expt.X)/expt.X) < 0.1)
    pch <- rep(pch.open, length(exptdata$T))
    pch[in10] <- pch.filled
    points(exptdata$T, exptdata[, prop], pch = pch)
  }
  # Plot regression line at a single P
  P <- mean(exptdata$P)
  message("EOSplot: plotting line for P = ", P, " bar")
  xs <- seq(xlim[1], xlim[2], length.out = 200)
  calc.X <- EOScalc(coefficients, xs, P, ...)
  lines(xs, calc.X, lty = lty, col = col)
  # Make legend
  if(!is.null(fun.legend) & !add) {
    # 20161101: negate QBorn and V_s_var
    iQ <- names(coefficients) %in% c("QBorn", "V_s_var")
    coefficients[iQ] <- -coefficients[iQ]
    coeffs <- as.character(round(as.numeric(coefficients), 4))
    # So that positive ones appear with a plus sign
    ipos <- which(coeffs >= 0)
    coeffs[ipos] <- paste("+", coeffs[ipos], sep = "")
    # Make labels for the functions
    fun.lab <- as.expression(lapply(1:length(coeffs),
      function(x) {EOSlab(names(coefficients)[x],coeffs[x])} ))
    #fun.lab <- paste(names(coeffs),round(as.numeric(coeffs),4))
    legend(fun.legend, legend = fun.lab, pt.cex = 0.1)
  }
  return(invisible(list(xrange = range(exptdata$T[iexpt]), coefficients = coefficients)))
}

EOScoeffs <- function(species, property, P = 1) {
  # Get the HKF coefficients for species in the database
  iis <- info(info(species, "aq"))
  if(property == "Cp") {
    out <- as.numeric(iis[,c("c1", "c2", "omega")])
    names(out) <- c("(Intercept)", "invTTheta2", "TXBorn")
  } else if(property == "V") {
    iis <- iis[,c("a1", "a2", "a3", "a4", "omega")]
    # Calculate sigma and xi and convert to volumetric units: 1 J = 10 cm^3 bar
    sigma <- convert( iis$a1 + iis$a2 / (2600 + P), "cm3bar" )
    xi <- convert( iis$a3 + iis$a4 / (2600 + P), "cm3bar" )
    omega <- convert( iis$omega, "cm3bar" )
    # 20151126: We _don't_ put a negative sign on omega here;
    # now, the negative sign in the HKF EOS is with the variable (QBorn or V_s_var)
    out <- c(sigma, xi, omega)
    names(out) <- c("(Intercept)", "invTTheta", "QBorn")
  }
  return(out)
}