Nothing
# 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.