[CHNOSZ-commits] r29 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 20 02:05:46 CET 2012
Author: jedick
Date: 2012-11-20 02:05:45 +0100 (Tue, 20 Nov 2012)
New Revision: 29
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
Log:
add invPPsi and invPPsiTTheta to EOSvar()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-11-11 11:46:47 UTC (rev 28)
+++ pkg/CHNOSZ/DESCRIPTION 2012-11-20 01:05:45 UTC (rev 29)
@@ -1,6 +1,6 @@
-Date: 2012-11-11
+Date: 2012-11-20
Package: CHNOSZ
-Version: 0.9.8-3
+Version: 0.9.8-4
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2012-11-11 11:46:47 UTC (rev 28)
+++ pkg/CHNOSZ/R/EOSregress.R 2012-11-20 01:05:45 UTC (rev 29)
@@ -14,6 +14,8 @@
"invTTheta" = (T-thermo$opt$Theta)^-1, # 1/(T-Theta)
"TTheta2" = (T-thermo$opt$Theta)^2, # (T-Theta)^2
"invTTheta2" = (T-thermo$opt$Theta)^-2, # 1/(T-Theta)^2
+ "invPPsi" = (P+thermo$opt$Psi)^-1, # 1/(P-Psi)
+ "invPPsiTTheta" = (P+thermo$opt$Psi)^-1 * (T-thermo$opt$Theta)^-1, # 1/[(P-Psi)(T-Theta)]
"V" = water(var,T=T,P=P)[,1],
"E" = water(var,T=T,P=P)[,1],
"kT" = water(var,T=T,P=P)[,1],
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2012-11-11 11:46:47 UTC (rev 28)
+++ pkg/CHNOSZ/inst/NEWS 2012-11-20 01:05:45 UTC (rev 29)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9.8-3 (2012-11-11)
+CHANGES IN CHNOSZ 0.9.8-4 (2012-11-20)
--------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
@@ -172,6 +172,10 @@
add today() for returning today's date in the format used in SUPCRT
files.
+- EOSvar() has new variables invPPsi and invPPsiTTheta, used for
+ temperature- and pressure-dependent regressions in the revised HKF
+ equations of state.
+
DATA UPDATES:
- Add 148 liquid and 148 crystalline acyclic isoprenoids, polycyclic
Modified: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2012-11-11 11:46:47 UTC (rev 28)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2012-11-20 01:05:45 UTC (rev 29)
@@ -44,6 +44,8 @@
\code{invTTheta} \tab \eqn{1/(T-\Theta)}{1/(T-Theta)} \cr
\code{TTheta2} \tab \eqn{(T-\Theta)^2}{(T-Theta)^2} \cr
\code{invTTheta2} \tab \eqn{1/(T-\Theta)^2}{1/(T-Theta)^2} \cr
+ \code{invPPsi} \tab \eqn{1/(P+\Psi)}{1/(P+Psi)} \cr
+ \code{invPPsiTTheta} \tab \eqn{1/((P+\Psi)(T-\Theta))}{1/((P+Psi)(T-Theta))} \cr
\code{V} \tab \eqn{V}{V} (volume of water) \cr
\code{E} \tab \eqn{E}{E} (isobaric expansivity of water) \cr
\code{kT} \tab \eqn{\kappa_T}{kT} (isothermal compressibility of water) \cr
@@ -84,61 +86,94 @@
\examples{
- \dontshow{data(thermo)}
- ## regress experimental heat capacities of CH4
- ## using revised Helgeson-Kirkham-Flowers equations
- # read the data from Hnedkovsky and Wood, 1997
- f <- system.file("extdata/cpetc/Cp.CH4.HW97.csv", package="CHNOSZ")
- d <- read.csv(f)
- # have to convert J to cal and MPa to bar
- d$Cp <- convert(d$Cp, "cal")
- d$P <- convert(d$P, "bar")
- # specify the terms in the HKF equations
- var <- c("invTTheta2", "TX")
- # perform regression, with a temperature limit
- EOSlm <- EOSregress(d, var, T.max=600)
- # the result is within 10% of the accepted
- # values of c1, c2 and omega for CH4(aq)
- CH4coeffs <- EOScoeffs("CH4", "Cp")
- dcoeffs <- EOSlm$coefficients - CH4coeffs
- stopifnot(all(abs(dcoeffs/CH4coeffs) < 0.1))
- ## make plots comparing the regressions
- ## here with the accepted EOS parameters of CH4
- par(mfrow=c(2,2))
- EOSplot(d, T.max=600)
- title("Cp of CH4(aq), fit to 600 K")
- legend("bottomleft", pch=1, legend="Hnedkovsky and Wood, 1997")
- EOSplot(d, coefficients=CH4coeffs)
- title("Cp from EOS parameters in database")
- EOSplot(d, T.max=600, T.plot=600)
- title("Cp fit to 600 K, plot to 600 K")
- EOSplot(d, coefficients=CH4coeffs, T.plot=600)
- title("Cp from EOS parameters in database")
+\dontshow{data(thermo)}
+## regress calculated heat capacities and volumes of CH4(aq)
+## to test that regressions return known HKF parameters
+# calculate the properties of CH4 on a T-P grid
+T <- convert(seq(0, 350, 25), "K")
+P <- seq(200, 1000, 100)
+# convert=FALSE means that temperature has units of K
+CH4.prop <- subcrt("CH4", T=T, P=P, grid="T", convert=FALSE)$out[[1]]
+# terms in the HKF equations for Cp
+Cp.var <- c("invTTheta2", "TX")
+# get coefficients in Cp regression
+Cp.lm <- EOSregress(CH4.prop[, c("T", "P", "Cp")], Cp.var)
+Cp.coeff <- Cp.lm$coefficients
+# terms in the HKF equations for V
+V.var <- c("invPPsi", "invTTheta", "invPPsiTTheta", "Q")
+# get coefficients in V regression
+V.lm <- EOSregress(CH4.prop[, c("T", "P", "V")], V.var)
+# use same units as HKF: convert from cm3.bar to calories (divide by 41.84)
+V.coeff <- convert(V.lm$coefficients, "calories")
+## the tests: did we get the HKF parameters that are in the database?
+CH4.par <- info(info("CH4"))
+# c1 and c2
+stopifnot(all.equal(Cp.coeff[1], CH4.par$c1, check.attr=FALSE))
+stopifnot(all.equal(Cp.coeff[2], CH4.par$c2, check.attr=FALSE))
+# omega (from Cp)
+stopifnot(all.equal(Cp.coeff[3], CH4.par$omega, check.attr=FALSE))
+# a1, a2, a3 and a4
+stopifnot(all.equal(V.coeff[1], CH4.par$a1, check.attr=FALSE))
+stopifnot(all.equal(V.coeff[2], CH4.par$a2, check.attr=FALSE))
+stopifnot(all.equal(V.coeff[3], CH4.par$a3, check.attr=FALSE))
+stopifnot(all.equal(V.coeff[4], CH4.par$a4, check.attr=FALSE))
+# omega (from V) - note negative sign
+stopifnot(all.equal(-V.coeff[5], CH4.par$omega, check.attr=FALSE))
- ## model experimental volumes of CH4
- ## using HKF equation and an exploratory one
- f <- system.file("extdata/cpetc/V.CH4.HWM96.csv", package="CHNOSZ")
- d <- read.csv(f)
- d$P <- convert(d$P, "bar")
- # the HKF equation
- varHKF <- c("invTTheta", "Q")
- # alpha is the expansivity coefficient of water
- varal <- c("invTTheta", "alpha")
- par(mfrow=c(2,2))
- # for both HKF and the expansivity equation
- # we'll fit up to a temperature limit
- EOSplot(d, varHKF, T.max=663, T.plot=625)
- legend("bottomright", pch=1, legend="Hnedkovsky et al., 1996")
- title("V of CH4(aq), HKF equation")
- EOSplot(d, varal, T.max=663, T.plot=625)
- title("V of CH4(aq), expansivity equation")
- EOSplot(d, varHKF, T.max=663)
- title("V of CH4(aq), HKF equation")
- EOSplot(d, varal, T.max=663)
- title("V of CH4(aq), expansivity equation")
- # note that the volume regression using the HKF gives
- # a result for omega (coefficient on Q) that is
- # not consistent with the high-T heat capacities
+## regress experimental heat capacities of CH4
+## using revised Helgeson-Kirkham-Flowers equations
+# read the data from Hnedkovsky and Wood, 1997
+f <- system.file("extdata/cpetc/Cp.CH4.HW97.csv", package="CHNOSZ")
+d <- read.csv(f)
+# have to convert J to cal and MPa to bar
+d$Cp <- convert(d$Cp, "cal")
+d$P <- convert(d$P, "bar")
+# specify the terms in the HKF equations
+var <- c("invTTheta2", "TX")
+# perform regression, with a temperature limit
+EOSlm <- EOSregress(d, var, T.max=600)
+# the result is within 10% of the accepted
+# values of c1, c2 and omega for CH4(aq)
+CH4coeffs <- EOScoeffs("CH4", "Cp")
+dcoeffs <- EOSlm$coefficients - CH4coeffs
+stopifnot(all(abs(dcoeffs/CH4coeffs) < 0.1))
+## make plots comparing the regressions
+## here with the accepted EOS parameters of CH4
+par(mfrow=c(2,2))
+EOSplot(d, T.max=600)
+title("Cp of CH4(aq), fit to 600 K")
+legend("bottomleft", pch=1, legend="Hnedkovsky and Wood, 1997")
+EOSplot(d, coefficients=CH4coeffs)
+title("Cp from EOS parameters in database")
+EOSplot(d, T.max=600, T.plot=600)
+title("Cp fit to 600 K, plot to 600 K")
+EOSplot(d, coefficients=CH4coeffs, T.plot=600)
+title("Cp from EOS parameters in database")
+
+## model experimental volumes of CH4
+## using HKF equation and an exploratory one
+f <- system.file("extdata/cpetc/V.CH4.HWM96.csv", package="CHNOSZ")
+d <- read.csv(f)
+d$P <- convert(d$P, "bar")
+# the HKF equation
+varHKF <- c("invTTheta", "Q")
+# alpha is the expansivity coefficient of water
+varal <- c("invTTheta", "alpha")
+par(mfrow=c(2,2))
+# for both HKF and the expansivity equation
+# we'll fit up to a temperature limit
+EOSplot(d, varHKF, T.max=663, T.plot=625)
+legend("bottomright", pch=1, legend="Hnedkovsky et al., 1996")
+title("V of CH4(aq), HKF equation")
+EOSplot(d, varal, T.max=663, T.plot=625)
+title("V of CH4(aq), expansivity equation")
+EOSplot(d, varHKF, T.max=663)
+title("V of CH4(aq), HKF equation")
+EOSplot(d, varal, T.max=663)
+title("V of CH4(aq), expansivity equation")
+# note that the volume regression using the HKF gives
+# a result for omega (coefficient on Q) that is
+# not consistent with the high-T heat capacities
}
\keyword{extra}
More information about the CHNOSZ-commits
mailing list