[CHNOSZ-commits] r104 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 27 15:15:54 CET 2015
Author: jedick
Date: 2015-11-27 15:15:54 +0100 (Fri, 27 Nov 2015)
New Revision: 104
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/demo/dehydration.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
pkg/CHNOSZ/tests/testthat/test-EOSregress.R
Log:
add ... argument to EOSvar() and EOSregress()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/DESCRIPTION 2015-11-27 14:15:54 UTC (rev 104)
@@ -1,6 +1,6 @@
-Date: 2015-11-25
+Date: 2015-11-27
Package: CHNOSZ
-Version: 1.0.7-1
+Version: 1.0.7-2
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/R/EOSregress.R 2015-11-27 14:15:54 UTC (rev 104)
@@ -3,7 +3,7 @@
# 20091105 first version
# 20110429 revise and merge with CHNOSZ package
-EOSvar <- function(var, T, P) {
+EOSvar <- function(var, T, P, ...) {
# get the variables of a term in a regression equation
# T (K), P (bar)
opt <- get("thermo")$opt
@@ -25,14 +25,17 @@
( if(var %in% water.props()) water(var, T, P)[, 1]
else if(exists(var)) {
if(is.function(get(var))) {
- if(identical(names(formals(get(var))), c("T", "P"))) get(var)(T, P)
- else stop(paste("the arguments of ", var, "() are not T, P", sep=""))
+ 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)
}
@@ -79,21 +82,21 @@
return(lab)
}
-EOSregress <- function(exptdata,var="",T.max=9999) {
+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,]
+ 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]
+ X <- exptdata[, 3]
# now build a regression formula
if(length(var) == 0) stop("var is missing")
- fmla <- as.formula(paste("X ~ ",paste(var,collapse="+")))
+ 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))
+ for(i in seq_along(var)) assign(var[i], EOSvar(var[i], T=T, P=P, ...))
# now regress away!
EOSlm <- lm(fmla)
return(EOSlm)
@@ -169,20 +172,21 @@
return(invisible(list(xlim=range(exptdata$T[iexpt]))))
}
-EOScoeffs <- function(species, property) {
+EOScoeffs <- function(species, property, P=1) {
# get the HKF coefficients for species in the database
iis <- info(info(species, "aq"))
if(property=="Cp") {
- out <- iis[,c("c1", "c2", "omega")]
+ 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 cal = 41.84 cm^3 bar
- sigma <- convert( iis$a1 + iis$a2 / (2600 + 1), "cm3bar" )
- xi <- convert( iis$a3 + iis$a4 / (2600 + 1), "cm3bar" )
+ sigma <- convert( iis$a1 + iis$a2 / (2600 + P), "cm3bar" )
+ xi <- convert( iis$a3 + iis$a4 / (2600 + P), "cm3bar" )
omega <- convert( iis$omega, "cm3bar" )
- # watch for the negative sign on omega here!
- out <- data.frame(sigma, xi, -omega)
+ # 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)
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/R/hkf.R 2015-11-27 14:15:54 UTC (rev 104)
@@ -17,7 +17,7 @@
property <- eargs$prop
props <- eargs$props
Prop <- eargs$Prop
- domega <- rep(domega,length.out=nrow(ghs))
+ domega <- rep(domega,length.out=nrow(eos))
# nonsolvation, solvation, and origination contribution
contribs <- c('n','s','o')
notcontrib <- ! contrib %in% contribs
@@ -29,8 +29,8 @@
# only take these ones if we're in SUPCRT92 compatibility mode
dosupcrt <- thermo$opt$water != "IAPWS95"
if(dosupcrt) {
- # (E, daldT, V - for partial derivatives of omega (g function))
- H2O.props <- c(H2O.props,'E','daldT','kT','ZBorn')
+ # (rho, alpha, daldT, beta - for partial derivatives of omega (g function))
+ H2O.props <- c(H2O.props, "rho", "alpha", "daldT", "beta")
} else {
# (NBorn, UBorn - for compressibility, expansibility)
H2O.props <- c(H2O.props,'NBorn','UBorn')
@@ -42,20 +42,26 @@
}
# a list to store the result
x <- list()
- for(k in 1:nrow(ghs)) {
+ nspecies <- nrow(ghs)
+ # we can be called with NULL ghs (by Cp_s_var, V_s_var)
+ if(is.null(nspecies)) nspecies <- nrow(eos)
+ for(k in 1:nspecies) {
# loop over each species
GHS <- ghs[k,]
EOS <- eos[k,]
# substitute Cp and V for missing EoS parameters
# here we assume that the parameters are in the same position as in thermo$obigt
- # put the heat capacity in for c1 if both c1 and c2 are missing
- if(all(is.na(EOS[, 17:18]))) EOS[, 17] <- EOS$Cp
- # put the volume in for a1 if a1, a2, a3 and a4 are missing
- if(all(is.na(EOS[, 13:16]))) EOS[, 13] <- convert(EOS$V, "calories")
- # test for availability of the EoS parameters
- hasEOS <- any(!is.na(EOS[, 13:20]))
- # if at least one of the EoS parameters is available, zero out any NA's in the rest
- if(hasEOS) EOS[, 13:20][, is.na(EOS[, 13:20])] <- 0
+ # we don't need this if we're just looking at solvation properties (Cp_s_var, V_s_var)
+ if("n" %in% contrib) {
+ # put the heat capacity in for c1 if both c1 and c2 are missing
+ if(all(is.na(EOS[, 17:18]))) EOS[, 17] <- EOS$Cp
+ # put the volume in for a1 if a1, a2, a3 and a4 are missing
+ if(all(is.na(EOS[, 13:16]))) EOS[, 13] <- convert(EOS$V, "calories")
+ # test for availability of the EoS parameters
+ hasEOS <- any(!is.na(EOS[, 13:20]))
+ # if at least one of the EoS parameters is available, zero out any NA's in the rest
+ if(hasEOS) EOS[, 13:20][, is.na(EOS[, 13:20])] <- 0
+ }
# compute values of omega(P,T) from those of omega(Pr,Tr)
# using g function etc. (Shock et al., 1992 and others)
omega <- EOS$omega # omega.PrTr
Modified: pkg/CHNOSZ/demo/dehydration.R
===================================================================
--- pkg/CHNOSZ/demo/dehydration.R 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/demo/dehydration.R 2015-11-27 14:15:54 UTC (rev 104)
@@ -1,6 +1,6 @@
# plot temperature dependence of log K for some dehydration reactions
-# the RSVGTipsDevice packages allows us to create an SVG file with
+# the RSVGTipsDevice package allows us to create an SVG file with
# tooltips and hyperlinks
if(require("RSVGTipsDevice")) {
@@ -14,7 +14,7 @@
plot(range(T), c(-2, 1), type="n", xlab="T, °C", ylab="log K")
title(main="Dehydration reactions")
-# add.obigt is needed to add malate and fumarate,
+# add.obigt is used to add malate and fumarate,
# and epsomite and hexahydrite to thermo$obigt
add.obigt()
reactants <- c("[AABB]", "[AABB]", "malate-2", "goethite", "gypsum", "epsomite", "ethanol")
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/inst/NEWS 2015-11-27 14:15:54 UTC (rev 104)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.7-1 (2015-11-25)
+CHANGES IN CHNOSZ 1.0.7-2 (2015-11-27)
--------------------------------------
- Add gypsum to OBIGT.csv.
@@ -9,7 +9,14 @@
- Add demo dehydration.R. This demo requires the RSVGTipsDevice package,
which is currently unavailable on Windows.
+- Add ... argument to EOSvar() and EOSregress(). Allows specifying new
+ variables that are written as a function of arbitrary properties
+ (regression variables still must be a function of at least T and P).
+- EOSvar() now returns the negative of the 'QBorn' variable so that
+ values of omega obtained from regressions of volume data have the
+ correct sign in the HKF equations.
+
CHANGES IN CHNOSZ 1.0.7 (2015-11-19)
------------------------------------
Modified: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2015-11-27 14:15:54 UTC (rev 104)
@@ -13,13 +13,13 @@
}
\usage{
- EOSregress(exptdata, var = "", T.max = 9999)
- EOSvar(var, T, P)
+ EOSregress(exptdata, var = "", T.max = 9999, ...)
+ EOSvar(var, T, P, ...)
EOScalc(coefficients, T, P)
EOSplot(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
fun.legend = "topleft", coefficients = NULL)
EOSlab(var, coeff = "")
- EOScoeffs(species, property)
+ EOScoeffs(species, property, P=1)
}
\arguments{
@@ -28,6 +28,7 @@
\item{T.max}{numeric, maximum temperature for regression, in degrees Kelvin}
\item{T}{numeric, temperature in degrees Kelvin}
\item{P}{numeric, pressure in bars}
+ \item{...}{arguments defining additional properties which variables are a function of}
\item{T.plot}{numeric, upper limit of temperature range to plot}
\item{fun.legend}{character, where to place legend on plot}
\item{coefficients}{dataframe, coefficients to use to make line on plot}
@@ -74,6 +75,7 @@
\code{EOScoeffs} retrieves coefficients in the Helgeson-Kirkham-Flowers equations from the thermodynamic database (\code{\link{thermo}$obigt}) for the given aqueous \code{species}.
If the \code{property} is \samp{Cp}, the resulting data frame has column names of \samp{(Intercept)}, \samp{invTTheta2} and \samp{TX}, respectively holding the coefficients \eqn{c_1}{c1}, \eqn{c_2}{c2} and \eqn{\omega}{omega} in the equation \eqn{Cp^\circ = c_1 + c_2/(T-\Theta)^2 + {\omega}TX}{Cp = c1 + c2/(T-Theta)^2 + omega*TX}.
If the \code{property} is \samp{V}, the data frame has column names of \samp{(Intercept)}, \samp{invTTheta} and \samp{Q}, respectively holding the coefficients \eqn{\sigma}{sigma}, \eqn{\xi}{xi} and \eqn{-\omega}{-omega} in \eqn{V^\circ = \sigma + \xi/(T-\Theta) - {\omega}Q}{V = sigma + xi/(T-Theta) - omega*Q}.
+Here, \eqn{\sigma}{sigma} and \eqn{\xi}{xi} are calculated from a1, a2, a3 and a4 in \code{thermo$obigt} at the pressure indicated by \code{P} (default 1 bar).
The motivation for writing these functions is to explore alternatives or possible modifications to the revised Helgeson-Kirkham-Flowers equations applied to aqueous nonelectrolytes.
As pointed out by Schulte et al., 2001, the functional forms of the equations do not permit retrieving values of the solvation parameter (\eqn{\omega}{omega}) that closely represent the observed trends in both heat capacity and volume at high temperatures (above ca. 200 \eqn{^{\circ}}{°}C).
Modified: pkg/CHNOSZ/tests/testthat/test-EOSregress.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-EOSregress.R 2015-11-25 04:43:03 UTC (rev 103)
+++ pkg/CHNOSZ/tests/testthat/test-EOSregress.R 2015-11-27 14:15:54 UTC (rev 104)
@@ -13,7 +13,7 @@
#expect_error(EOSvar("TX", T=25, P=1), "the arguments of TX\\(\\) are not T, P")
})
-test_that("regressions return known HKF parameters", {
+test_that("regressions return known HKF parameters (neutral species)", {
# regress computed values of heat capacity and volume of CH4(aq)
# calculated from HKF parameters on a T-P grid
T <- convert(seq(0, 350, 25), "K")
@@ -34,15 +34,15 @@
## the tests: did we get the HKF parameters that are in the database?
CH4.par <- info(info("CH4"))
# c1 and c2
- expect_equal(Cp.coeff[1], CH4.par$c1, check.attributes=FALSE)
- expect_equal(Cp.coeff[2], CH4.par$c2, check.attributes=FALSE)
+ expect_equivalent(Cp.coeff[1], CH4.par$c1)
+ expect_equivalent(Cp.coeff[2], CH4.par$c2)
# omega (from Cp)
- expect_equal(Cp.coeff[3], CH4.par$omega, check.attributes=FALSE)
+ expect_equivalent(Cp.coeff[3], CH4.par$omega)
# a1, a2, a3 and a4
- expect_equal(V.coeff[1], CH4.par$a1, check.attributes=FALSE)
- expect_equal(V.coeff[2], CH4.par$a2, check.attributes=FALSE)
- expect_equal(V.coeff[3], CH4.par$a3, check.attributes=FALSE)
- expect_equal(V.coeff[4], CH4.par$a4, check.attributes=FALSE)
- # omega (from V) - note negative sign
- expect_equal(-V.coeff[5], CH4.par$omega, check.attributes=FALSE)
+ expect_equivalent(V.coeff[1], CH4.par$a1)
+ expect_equivalent(V.coeff[2], CH4.par$a2)
+ expect_equivalent(V.coeff[3], CH4.par$a3)
+ expect_equivalent(V.coeff[4], CH4.par$a4)
+ # omega (from V)
+ expect_equivalent(V.coeff[5], CH4.par$omega)
})
More information about the CHNOSZ-commits
mailing list