[CHNOSZ-commits] r42 - in pkg/CHNOSZ: . R data demo inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 12 23:35:25 CET 2013
Author: jedick
Date: 2013-02-12 23:35:25 +0100 (Tue, 12 Feb 2013)
New Revision: 42
Added:
pkg/CHNOSZ/inst/tests/test-EOSregress.R
pkg/CHNOSZ/inst/tests/test-water.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/CO2Ac.R
pkg/CHNOSZ/demo/nonideal.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/water.Rd
Log:
use consistent names for properties of water
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/DESCRIPTION 2013-02-12 22:35:25 UTC (rev 42)
@@ -1,6 +1,6 @@
-Date: 2013-02-06
+Date: 2013-02-13
Package: CHNOSZ
-Version: 0.9-9.4
+Version: 0.9-9.5
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 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/R/EOSregress.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -3,11 +3,11 @@
# 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)
out <- switch(EXPR = var,
- "(Intercept)" = rep(1,length(T)),
+ "(Intercept)" = rep(1, length(T)),
"T" = T,
"P" = P,
"TTheta" = T-thermo$opt$Theta, # T-Theta
@@ -16,17 +16,12 @@
"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],
- "alpha" = water(var,T=T,P=P)[,1],
- "beta" = water(var,T=T,P=P)[,1],
- "X" = water(var,T=T,P=P)[,1],
- "Q" = water(var,T=T,P=P)[,1],
- "TX" = T*water("X",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],
- NA
+ "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],
+ # the "default": get a variable that is a property of water
+ (if(var %in% water.props()) water(var, T, P)[, 1]
+ else stop(paste("can't find a variable named", var)))
)
return(out)
}
@@ -43,9 +38,9 @@
"P" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
"V" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
"E" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "X" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "Q" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "TX" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
+ "XBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
+ "QBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
+ "TXBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
"kT" = substitute(YYY%*%kappa[italic(T)],list(YYY=coeff)),
"alpha" = substitute(YYY%*%alpha,list(YYY=coeff)),
"beta" = substitute(YYY%*%beta,list(YYY=coeff)),
@@ -95,8 +90,8 @@
prop <- colnames(exptdata)[3]
# if var is NULL use HKF equations
if(is.null(var)) {
- if(prop=="Cp") var <- c("invTTheta2","TX")
- if(prop=="V") var <- c("invTTheta","Q")
+ if(prop=="Cp") var <- c("invTTheta2","TXBorn")
+ if(prop=="V") var <- c("invTTheta","QBorn")
}
expt <- exptdata
# perform the regression, only using temperatures up to T.max
@@ -165,19 +160,19 @@
return(invisible(list(xlim=range(expt$T[iexpt]))))
}
-EOScoeffs <- function(species,property) {
+EOScoeffs <- function(species, property) {
# get the HKF coefficients for species in the database
- iis <- info(info(species,"aq"))
+ iis <- info(info(species, "aq"))
if(property=="Cp") {
- out <- iis[,c("c1","c2","omega")]
- names(out) <- c("(Intercept)","invTTheta2","TX")
+ out <- iis[,c("c1", "c2", "omega")]
+ names(out) <- c("(Intercept)", "invTTheta2", "TXBorn")
} else if(property=="V") {
- iis <- iis[,c("a1","a2","a3","a4","omega")]
+ iis <- iis[,c("a1", "a2", "a3", "a4", "omega")]
sigma <- ( iis$a1 + iis$a2 / (2600 + 1) ) * 41.84
xi <- ( iis$a3 + iis$a4 / (2600 + 1) ) * 41.84
# watch for the negative sign on omega here!
- out <- data.frame(sigma,xi,-iis$omega)
- names(out) <- c("(Intercept)","invTTheta","Q")
+ out <- data.frame(sigma, xi, -iis$omega)
+ names(out) <- c("(Intercept)", "invTTheta", "QBorn")
}
return(out)
}
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/R/hkf.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -24,20 +24,20 @@
stop(paste('argument',c2s(contrib[notcontrib]),'not in',c2s(contribs),'n'))
# get water properties, if they weren't supplied in arguments (and we want solvation props)
if('s' %in% contrib) {
- H2O.props <- c('Q','X','Y','epsilon')
+ H2O.props <- c("QBorn", "XBorn", "YBorn", "diel")
# only take these ones if we're in SUPCRT92 compatibility mode
- dosupcrt <- length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3))!=0
+ 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','Z')
+ H2O.props <- c(H2O.props,'E','daldT','kT','ZBorn')
} else {
- # (N, UBorn - for compressibility, expansibility)
- H2O.props <- c(H2O.props,'N','UBorn')
+ # (NBorn, UBorn - for compressibility, expansibility)
+ H2O.props <- c(H2O.props,'NBorn','UBorn')
}
if(is.null(H2O.PT)) H2O.PT <- water(H2O.props,T=T,P=P)
if(is.null(H2O.PrTr)) H2O.PrTr <- water(H2O.props,T=thermo$opt$Tr,P=thermo$opt$Pr)
- ZBorn <- -1/H2O.PT$epsilon
- ZBorn.PrTr <- -1/H2O.PrTr$epsilon
+ ZBorn <- -1/H2O.PT$diel
+ ZBorn.PrTr <- -1/H2O.PrTr$diel
}
# a list to store the result
x <- list()
@@ -125,20 +125,20 @@
if( icontrib=="s") {
# solvation ghs equations
if(prop=="g") {
- p <- -omega.PT*(ZBorn+1) + omega*(ZBorn.PrTr+1) + omega*H2O.PrTr$Y*(T-Tr)
+ p <- -omega.PT*(ZBorn+1) + omega*(ZBorn.PrTr+1) + omega*H2O.PrTr$YBorn*(T-Tr)
# at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA
if(!is.na(GHS$G)) p[T==Tr & P==Pr] <- 0
}
if(prop=="h")
- p <- -omega.PT*(ZBorn+1) + omega.PT*T*H2O.PT$Y + T*(ZBorn+1)*dwdT +
- omega*(ZBorn.PrTr+1) - omega*Tr*H2O.PrTr$Y
+ p <- -omega.PT*(ZBorn+1) + omega.PT*T*H2O.PT$YBorn + T*(ZBorn+1)*dwdT +
+ omega*(ZBorn.PrTr+1) - omega*Tr*H2O.PrTr$YBorn
if(prop=="s")
- p <- omega.PT*H2O.PT$Y + (ZBorn+1)*dwdT - omega*H2O.PrTr$Y
+ p <- omega.PT*H2O.PT$YBorn + (ZBorn+1)*dwdT - omega*H2O.PrTr$YBorn
# solvation cp v kt e equations
- if(prop=='cp') p <- omega.PT*T*H2O.PT$X + 2*T*H2O.PT$Y*dwdT +
+ if(prop=='cp') p <- omega.PT*T*H2O.PT$XBorn + 2*T*H2O.PT$YBorn*dwdT +
T*(ZBorn+1)*d2wdT2
if(prop=='v') p <- -convert(omega.PT,'cm3bar') *
- H2O.PT$Q + convert(dwdP,'cm3bar') * (-ZBorn - 1)
+ H2O.PT$QBorn + convert(dwdP,'cm3bar') * (-ZBorn - 1)
# WARNING: the partial derivatives of omega are not included here here for kt and e
# (to do it, see p. 820 of SOJ+92 ... but kt requires d2wdP2 which we don't have yet)
if(prop=='kt') p <- convert(omega,'cm3bar') * H2O.PT$N
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/R/subcrt.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -267,16 +267,16 @@
# than possible many times in hkf()).
wprop.PT <- character()
wprop.PrTr <- 'rho'
- dosupcrt <- length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3))!=0
- if(TRUE %in% (prop %in% c('logk','g','h','s'))) wprop.PrTr <- c(wprop.PrTr,'Y')
- if(dosupcrt | TRUE %in% (prop %in% c('logk','g','h'))) wprop.PrTr <- c(wprop.PrTr,'epsilon')
+ dosupcrt <- thermo$opt$water != "IAPWS95"
+ if(TRUE %in% (prop %in% c('logk','g','h','s'))) wprop.PrTr <- c(wprop.PrTr,'YBorn')
+ if(dosupcrt | TRUE %in% (prop %in% c('logk','g','h'))) wprop.PrTr <- c(wprop.PrTr,'diel')
H2O.PrTr <- water(wprop.PrTr,T=thermo$opt$Tr,P=thermo$opt$Pr)
- if(TRUE %in% (prop %in% c('cp'))) {wprop.PT <- c(wprop.PT,'X','Y')}
- if(TRUE %in% (prop %in% c('v'))) {wprop.PT <- c(wprop.PT,'Q')}
- if(TRUE %in% (prop %in% c('kt'))) {wprop.PT <- c(wprop.PT,'N')}
+ if(TRUE %in% (prop %in% c('cp'))) {wprop.PT <- c(wprop.PT,'XBorn','YBorn')}
+ if(TRUE %in% (prop %in% c('v'))) {wprop.PT <- c(wprop.PT,'QBorn')}
+ if(TRUE %in% (prop %in% c('kt'))) {wprop.PT <- c(wprop.PT,'NBorn')}
if(TRUE %in% (prop %in% c('e'))) {wprop.PT <- c(wprop.PT,'UBorn')}
# get additional properties required for omega derivatives
- if(dosupcrt) wprop.PT <- c(wprop.PT,'alpha','daldT','beta','epsilon')
+ if(dosupcrt) wprop.PT <- c(wprop.PT,'alpha','daldT','beta','diel')
H2O.PT <- water(c(wprop.PrTr,wprop.PT),T=T,P=P)
if(TRUE %in% isaq) {
# now the species stuff
Modified: pkg/CHNOSZ/R/util.args.R
===================================================================
--- pkg/CHNOSZ/R/util.args.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/R/util.args.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -31,23 +31,21 @@
return(list(props=props,prop=prop,Prop=Prop))
}
-TP.args <- function(T=NULL,P=NULL) {
- if(!is.null(P)) {
- if(identical(P[1],'Psat')) {
- P <- water('Psat',T,P=NULL)
- P <- P[,1]
- # water.SUPCRT92 issues its own warnings about
- # exceeding Psat's temperature limit
- if(length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3))==0)
- if(length(which(is.na(P)))>0)
- warning('TP.args: NAs in Psat (likely T > Tc where Tc = 647.096 K)',call.=FALSE)
- }
+TP.args <- function(T=NULL, P=NULL) {
+ # keep the [1] here because some functions (e.g. subcrt) will repeat "Psat"
+ if(identical(P[1], "Psat")) {
+ P <- water("Psat", T, P="Psat")[, 1]
+ # water.SUPCRT92 issues its own warnings about
+ # exceeding Psat's temperature limit
+ if(thermo$opt$water == "IAPWS95")
+ if(length(which(is.na(P)))>0)
+ warning('TP.args: NAs in Psat (likely T > Tc where Tc = 647.096 K)',call.=FALSE)
}
if(length(P) < length(T) & !is.null(P)) P <- rep(P, length.out=length(T))
else if(length(T) < length(P) & !is.null(T)) T <- rep(T, length.out=length(P))
# something we do here so the SUPCRT water calculations work
T[T==273.15] <- 273.16
- return(list(T=T,P=P))
+ return(list(T=T, P=P))
}
state.args <- function(state=NULL) {
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/R/water.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -51,207 +51,183 @@
return(t)
}
-water <- function(property = NULL,T = thermo$opt$Tr, P = 'Psat') {
+water <- function(property = NULL, T = thermo$opt$Tr, P = "Psat") {
# calculate the properties of liquid H2O as a function of T and P
# T in Kelvin, P in bar
if(is.null(property)) stop('property was NULL')
- # this tells us to do the calculations using code taken from SUPCRT
- do.supcrt <- length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3)) > 0
- eargs <- eos.args('water',property=property,T=T,P=P)
- property <- eargs$prop; Property <- eargs$Prop
- # working out the arguments
- tpargs <- TP.args(T=T,P=P)
- P <- tpargs$P
- T <- tpargs$T
- # Psat stuff
- psat <- function() {
- p <- numeric()
- if(do.supcrt) {
- p <- water.SUPCRT92('',T=T,rep(0,length(T)),isat=1)
- p[p==0] <- NaN
- return(p)
- } else {
- for(i in 1:length(T)) {
- if(T[i] < 373.124) p <- c(p,0.1)
- else p <- c(p,WP02.auxiliary('P.sigma',T[i]))
- }
- return(convert(p,'bar'))
- }
+ # make T and P equal length
+ if(!identical(P, "Psat")) {
+ if(length(P) < length(T)) P <- rep(P, length.out = length(T))
+ else if(length(T) < length(P)) T <- rep(T, length.out = length(P))
}
- # a quick return if property = 'Psat', for use by the TP.args() function
- if(length(property)==1 & property[1]=='psat') return(data.frame(Psat=psat()))
- ### maybe we are using the SUPCRT calculations ###
+ # turn 273.15 K to 273.16 K (needed for water.SUPCRT92 at Psat)
+ T[T == 273.15] <- 273.16
+ # this tells us to do the calculations using code taken from SUPCRT
+ do.supcrt <- thermo$opt$water != "IAPWS95"
if(do.supcrt) {
- names.SUPCRT <- c('Speed','alpha','beta','alpha','beta','diel','ZBorn','YBorn','QBorn','XBorn')
- names.CHNOSZ <- c('w','alpha','beta','E','kT','epsilon','Z','Y','Q','X')
- Property.new <- character()
- # convert names to SUPCRT
- for(i in 1:length(Property)) if(Property[i] %in% names.CHNOSZ)
- Property.new[i] <- names.SUPCRT[match(Property[i],names.CHNOSZ)]
- else Property.new[i] <- Property[i]
- # deal with compressibility and expansivity 20091203
- iE <- which(Property=="E")
- ikT <- which(Property=="kT")
- iV <- numeric()
- if("kT" %in% Property | "E" %in% Property) iV <- length(Property.new <- c(Property.new,"V"))
- # get the value of the property
- w.out <- water.SUPCRT92(Property.new,T=T,P=P)
- # finish dealing with compressibility and expansivity
- if("E" %in% Property) w.out[,iE] <- w.out$V*w.out$alpha
- if("kT" %in% Property) w.out[,ikT] <- w.out$V*w.out$beta
- if(length(iV) > 0) w.out <- w.out[,-iV,drop=FALSE]
- colnames(w.out) <- Property
+ # get the values of the properties using SUPCRT92
+ w.out <- water.SUPCRT92(property, T, P)
return(w.out)
} else {
# here we get properties using IAPWS-95
w.out <- water.IAPWS95(property, T, P)
- colnames(w.out) <- Property
+ colnames(w.out) <- property
return(w.out)
}
}
+water.props <- function(formulation=thermo$opt$water) {
+ # return the names of properties that are available in SUPCRT92 or IAPWS95
+ # added 20130212 jmd
+ if(formulation=="SUPCRT92")
+ props <- c("A", "G", "S", "U", "H", "Cv", "Cp",
+ "Speed", "alpha", "beta", "diel", "visc",
+ "tcond", "surten", "tdiff", "Prndtl", "visck", "albe",
+ "ZBorn", "YBorn", "QBorn", "daldT", "XBorn",
+ "V", "rho", "Psat", "E", "kT")
+ else if(formulation=="IAPWS95")
+ props <- c("A", "G", "S", "U", "H", "Cv", "Cp",
+ "YBorn", "QBorn", "XBorn", "NBorn", "UBorn",
+ "V", "rho", "Psat", "de.dT", "de.dP", "P")
+ return(props)
+}
-water.SUPCRT92 <- function(property,T=298.15,P=1,isat=0) {
+water.SUPCRT92 <- function(property, T=298.15, P=1) {
### interface to H2O92D.f : FORTRAN subroutine taken from
### SUPCRT92 for calculating the thermodynamic and
### electrostatic properties of H2O.
## we restrict the calculations to liquid water
## except for getting Psat (vapor-liquid saturation
## pressure as a function of T>100 C). 20071213 jmd
- # H2O92 doesn't output Born functions N or U
- if('n' %in% tolower(property) | 'uborn' %in% tolower(property))
- stop('I can\'t tell you the Born functions N or U (used in calculating compressibilities and expansibilities of aqueous species).')
- # pressure setting
- if(is.null(P)) P <- rep(0,length(T))
- # values to use here gleaned from H2O92D.f and SUP92D.f
+ # check for availability of properties
+ iprop <- match(tolower(property), tolower(water.props("SUPCRT92")))
+ if(any(is.na(iprop))) stop(paste("property(s) not available:", paste(property[is.na(iprop)], collapse=" ")))
+ # make sure Psat in properties comes with isat=1
+ if("psat" %in% tolower(property) & !identical(P, "Psat")) stop("please set P='Psat' to calculate the property Psat")
+ # for Psat(T) (1) or T-P (2)
+ if(identical(P, "Psat")) iopt <- 1 else iopt <- 2
+ if(identical(P, "Psat")) isat <- 1 else isat <- 0
+ # input values, gleaned from H2O92D.f and SUP92D.f
# it, id, ip, ih, itripl, isat, iopt, useLVS, epseqn, icrit
- if(isat) iopt <- 1 else iopt <- 2 # for Psat(T) (1) or T-P (2)
- specs <- c(2,2,2,5,1,isat,iopt,1,4,0)
- states <- rep(0,4)
- # match up properties with the output
- props <- c('a','g','s','u','h','cv','cp','Speed','alpha',
- 'beta','diel','visc','tcond','surten','tdiff','Prndtl',
- 'visck','albe','ZBorn','YBorn','QBorn','daldT','XBorn')
- iprop <- seq(1,45,length.out=23)
- # now to the actual calculations
- Tc <- convert(T,'C')
+ specs <- c(2, 2, 2, 5, 1, isat, iopt, 1, 4, 0)
+ states <- rep(0, 4)
# initialize the output matrix
- w.out <- matrix(NA,nrow=length(T),ncol=23,byrow=TRUE)
+ w.out <- matrix(NA, nrow=length(T), ncol=23, byrow=TRUE)
err.out <- numeric(length(T))
rho.out <- numeric(length(T))
- p.out <- numeric(length(T))
+ P.out <- numeric(length(T))
# 20091022 TODO: parallelize this
+ Tc <- convert(T, "C")
for(i in 1:length(T)) {
states[1] <- Tc[i]
- states[2] <- P[i]
- if(any(is.na(c(Tc[i],P[i])))) {
+ if(identical(P, "Psat")) states[2] <- 0
+ else states[2] <- P[i]
+ if(is.na(Tc[i]) | is.na(P[i]) & !identical(P, "Psat")) {
# if T or P is NA, all properties are NA
- w <- matrix(rep(NA,23),nrow=1)
- w.out[i,] <- w
- p.out[i] <- NA
+ # (NA's are already in w.out)
+ P.out[i] <- NA
rho.out[i] <- NA
} else {
- inc <- 0
- h2o <- .Fortran('H2O92',as.integer(specs),as.double(states),
- as.double(rep(0,46)),as.integer(0),PACKAGE='CHNOSZ')
+ # now to the actual calculations
+ H2O <- .Fortran("H2O92", as.integer(specs), as.double(states),
+ as.double(rep(0, 46)), as.integer(0), PACKAGE="CHNOSZ")
# errors
- err <- h2o[[4]]
- err.out[i] <- err
- # density
- rho <- h2o[[2]][3]
- rho2 <- h2o[[2]][4]
+ err.out[i] <- H2O[[4]]
+ # density of two states
+ rho <- H2O[[2]][3]
+ rho2 <- H2O[[2]][4]
if(rho2 > rho) {
# liquid is denser than vapor
rho <- rho2
- # for selecting the liquid properties later
- inc <- 1
- }
+ inc <- 1 # second state is liquid
+ } else inc <- 0 # first state is liquid
rho.out[i] <- rho
- # most of the properties we're interested in
- w <- t(h2o[[3]][iprop+inc])
- if(err==1) w[1,] <- NA
+ # 23 properties of the phase in the liquid state
+ w <- t(H2O[[3]][seq(1, 45, length.out=23)+inc])
+ if(err.out[i]==1) w[1, ] <- NA
# update the ith row of the output matrix
w.out[i,] <- w
# Psat
- if(isat | 'psat' %in% tolower(property)) {
- p <- h2o[[2]][2]
- p[p==0] <- NA
+ if(identical(P, "Psat")) {
+ w.P <- H2O[[2]][2]
+ w.P[w.P==0] <- NA
# Psat specifies P=1 below 100 degC
- p[p < 1] <- 1
- p.out[i] <- p
- } else {
- p.out[i] <- P[i]
+ w.P[w.P < 1] <- 1
+ P.out[i] <- w.P
}
}
}
# convert output to dataframe
w.out <- as.data.frame(w.out)
- names(w.out) <- props
- # assemble the properties
- mwH2O <- 18.0152 # SUP92.f
- w.out <- cbind(w.out,V=mwH2O/rho.out,rho=rho.out*1000)
- if(isat | 'psat' %in% tolower(property)) w.out <- cbind(w.out,Psat=p.out)
+ # add names of properties to the output
+ names(w.out) <- water.props("SUPCRT92")[1:23]
+ # assemble additional properties: V, rho, Psat, E, kT
+ if(any(iprop > 23)) {
+ mwH2O <- 18.0152 # SUP92.f
+ V=mwH2O/rho.out
+ rho=rho.out*1000
+ Psat=P.out
+ E <- V*w.out$alpha
+ kT <- V*w.out$beta
+ w.out <- cbind(w.out, data.frame(V=V, rho=rho, Psat=Psat, E=E, kT=kT))
+ }
# tell the user about any problems
if(any(err.out==1)) {
if(length(T) > 1) plural <- "s" else plural <- ""
nerr <- length(which(err.out==1))
if(nerr > 1) plural2 <- "s" else plural2 <- ""
- if(isat) msgout(paste("water.SUPCRT92: error",plural2," calculating ",
- nerr," of ",length(T)," point",plural,"; for Psat we need 273.16 < T < 647.067 K\n",sep=""))
- else msgout(paste("water.SUPCRT92: error",plural2," calculating ",nerr,
- " of ",length(T)," point",plural,
- "; T < Tfusion at P, T > 2250 degC, or P > 30kb.\n",sep=""))
- # that last bit is taken from SUP92D.f in the SUPCRT92 distribution
+ if(identical(P, "Psat")) msgout(paste("water.SUPCRT92: error", plural2, " calculating ",
+ nerr, " of ", length(T), " point", plural, "; for Psat we need 273.16 < T < 647.067 K\n", sep=""))
+ else msgout(paste("water.SUPCRT92: error", plural2, " calculating ", nerr,
+ " of ", length(T), " point", plural,
+ "; T < Tfusion at P, T > 2250 degC, or P > 30kb.\n", sep=""))
+ # that last bit is taken from SUP92D.f in SUPCRT92
}
- # if isat is 1, just return the calculated pressures
- if(isat) return(w.out$Psat)
# return only the selected properties
- icol <- match(tolower(property),tolower(colnames(w.out)))
- return(w.out[,icol,drop=FALSE])
+ return(w.out[, iprop, drop=FALSE])
}
-water.IAPWS95 <- function(property, T=298.15, P=1, quiet=FALSE) {
- # to get the properties of water via IAPWS-95
- if(!quiet) msgout(paste("water.IAPWS95: calculating", length(T), "values for"))
- M <- 18.015268 # g mol-1
- rho <- function() {
- # return a density in kg m-3
- # corresponding to the given pressure (MPa) and temperature (K)
- pfun <- function(rho,T,P) {
- P <- convert(P,'MPa')
- t <- IAPWS95('p',rho=rho,T=T)[,1] - P
- return(t)
+rho.IAPWS95 <- function(T=298.15, P=1) {
+ # return a density in kg m-3
+ # corresponding to the given pressure (bar) and temperature (K)
+ if(identical(P, "Psat")) stop("this function doesn't take P='Psat'")
+ dP <- function(rho, T, P) {
+ dP <- IAPWS95("P", rho=rho, T=T)[, 1] - convert(P, "MPa")
+ return(dP)
+ }
+ rho <- numeric()
+ for(i in 1:length(T)) {
+ if(T[i] < 647.096) {
+ rho.lower <- WP02.auxiliary("rho.liquid", T=T[i]) - 2
+ rho.upper <- rho.lower + 400
+ if(P[i] < 5000) rho.upper <- rho.lower + 300
+ if(P[i] < 1000) rho.upper <- rho.lower + 200
+ if(P[i] < 300) rho.upper <- rho.lower + 30
+ } else {
+ rho.lower <- 0.01
+ rho.upper <- 1200
}
- t <- numeric()
- for(i in 1:length(T)) {
- if(T[i] < 647.096) {
- rho.lower <- WP02.auxiliary('rho.liquid',T=T[i])-2
- rho.upper <- rho.lower + 400
- if(P[i] < 5000) rho.upper <- rho.lower + 300
- if(P[i] < 1000) rho.upper <- rho.lower + 200
- if(P[i] < 300) rho.upper <- rho.lower + 30
- }
- else { rho.lower <- 0.01; rho.upper <- 1200}
- tu <- try(uniroot(pfun,c(rho.lower,rho.upper),T=T[i],P=P[i])$root,TRUE)
- if(!is.numeric(tu)) {
- warning('water: no root for density between ',round(rho.lower,1),
- ' and ',round(rho.upper,1),' kg m-3 at ',T[i],' K and ',P[i],' bar.',call.=FALSE,immediate.=TRUE)
- tu <- NA
- }
- t <- c(t,tu)
+ this.rho <- try(uniroot(dP, c(rho.lower, rho.upper), T=T[i], P=P[i])$root, TRUE)
+ if(!is.numeric(this.rho)) {
+ warning("rho.IAPWS95: no root for density between ", round(rho.lower, 1),
+ " and ", round(rho.upper, 1), " kg m-3 at ", T[i], " K and ", P[i], " bar", call.=FALSE)
+ this.rho <- NA
}
- return(t)
+ rho <- c(rho, this.rho)
}
+ return(rho)
+}
+
+water.IAPWS95 <- function(property, T=298.15, P=1) {
+ # to get the properties of water via IAPWS-95
+ msgout(paste("water.IAPWS95: calculating", length(T), "values for"))
+ M <- 18.015268 # g mol-1
v <- function() return(M*1000/my.rho)
- p <- function() return(P)
# Psat stuff
psat <- function() {
- p <- numeric()
- for(i in 1:length(T)) {
- if(T[i] < 373.124) p <- c(p,0.1)
- else p <- c(p,WP02.auxiliary('P.sigma',T[i]))
- }
- return(convert(p,'bar'))
+ p <- WP02.auxiliary("P.sigma", T)
+ p[T < 373.124] <- 0.1
+ return(convert(p, "bar"))
}
## thermodynamic properties
# convert to SUPCRT reference state
@@ -263,6 +239,8 @@
# does the reference state used for GHS also go here?
dU <- -67434.5 - 451.3229
dA <- -55814.06 + 20.07376 - dS * (T - thermo$opt$Tr)
+ # calculate pressure from the given T and estimated rho
+ p <- function() return(convert(IAPWS95("p", T=T, rho=my.rho), "bar"))
# convert IAPWS95() (specific, joule) to (molar, cal)
s <- function()
return(convert(IAPWS95('s',T=T,rho=my.rho)$s*M,'cal')+dS)
@@ -292,7 +270,7 @@
this.P <- P[i]
this.rho <- my.rho[i]
dt <- 0.001; t1 <- this.T-dt; t2 <- this.T+dt
- rho <- water.IAPWS95("rho", T=c(t1, t2), P=this.P, quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=c(t1, t2), P=rep(this.P, 2))
e <- water.AW90(T=c(t1,t2),rho=rho,rep(this.P,2))
p <- c(p,(e[2]-e[1])/(2*dt))
}
@@ -305,54 +283,54 @@
this.P <- P[i]
this.rho <- my.rho[i]
dp <- 0.001; p1 <- this.P-dp; p2 <- this.P+dp
- rho <- water.IAPWS95("rho", T=this.T, P=c(p1, p2), quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=rep(this.T, 2), P=c(p1, p2))
e <- water.AW90(P=c(p1,p2),rho=rho,T=rep(this.T,2))
p <- c(p,(e[2]-e[1])/(2*dp))
}
return(p)
}
## Born functions
- q <- function() {
+ qborn <- function() {
p <- numeric()
for(i in 1:length(T)) {
this.T <- T[i]; this.P <- P[i]; this.rho <- my.rho[i]
dp <- 0.01; p1 <- this.P-dp; p2 <- this.P+dp
- rho <- water.IAPWS95("rho", T=rep(this.T, 2), P=c(p1, p2), quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=rep(this.T, 2), P=c(p1, p2))
e <- water.AW90(T=rep(this.T,2),rho=rho,P=convert(c(p1,p2),'MPa'))
#p <- c(p,convert(-(1/e[2]-1/e[1])/(2*dp),'cm3bar'))
p <- c(p,-(1/e[2]-1/e[1])/(2*dp))
}
return(p)
}
- n <- function() {
+ nborn <- function() {
p <- numeric()
for(i in 1:length(T)) {
this.T <- T[i]; this.P <- P[i]; this.rho <- my.rho[i]
dp <- 0.01; p1 <- this.P-dp; p2 <- this.P+dp
- rho <- water.IAPWS95("rho", T=rep(this.T, 3), P=c(p1, this.P, p2), quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=rep(this.T, 3), P=c(p1, this.P, p2))
e <- water.AW90(T=rep(this.T,3),rho=rho,P=convert(c(p1,this.P,p2),'MPa'))
#p <- c(p,convert(convert((-(1/e[3]-1/e[2])/dp+(1/e[2]-1/e[1])/dp)/dp,'cm3bar'),'cm3bar'))
p <- c(p,(-(1/e[3]-1/e[2])/dp+(1/e[2]-1/e[1])/dp)/dp)
}
return(p)
}
- y <- function() {
+ yborn <- function() {
p <- numeric()
for(i in 1:length(T)) {
this.T <- T[i]; this.P <- P[i]; this.rho <- my.rho[i]
dt <- 0.001; t1 <- this.T-dt; t2 <- this.T+dt
- rho <- water.IAPWS95("rho", T=c(t1, t2), P=rep(this.P, 2), quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=c(t1, t2), P=rep(this.P, 2))
e <- water.AW90(T=c(t1,t2),rho=rho,P=convert(rep(this.P,2),'MPa'))
p <- c(p,-(1/e[2]-1/e[1])/(2*dt))
}
return(p)
}
- x <- function() {
+ xborn <- function() {
p <- numeric()
for(i in 1:length(T)) {
this.T <- T[i]; this.P <- P[i]; this.rho <- my.rho[i]
dt <- 0.001; t1 <- this.T-dt; t2 <- this.T+dt
- rho <- water.IAPWS95("rho", T=c(t1, this.T, t2), P=rep(this.P, 3), quiet=TRUE)[, 1]
+ rho <- rho.IAPWS95(T=c(t1, this.T, t2), P=rep(this.P, 3))
e <- water.AW90(T=c(t1,this.T,t2),rho=rho,P=convert(rep(this.P,3),'MPa'))
p <- c(p,(-(1/e[3]-1/e[2])/dt+(1/e[2]-1/e[1])/dt)/dt)
}
@@ -364,8 +342,8 @@
this.T <- T[i]; this.P <- P[i]; this.rho <- my.rho[i]
dt <- 0.001; this.T1 <- this.T - dt; this.T2 <- this.T + dt
dp <- 0.001; p1 <- this.P-dp; p2 <- this.P+dp
- rho1 <- water.IAPWS95("rho", T=rep(this.T1, 2), P=c(p1, p2), quiet=TRUE)[, 1]
- rho2 <- water.IAPWS95("rho", T=rep(this.T2, 2), P=c(p1, p2), quiet=TRUE)[, 1]
+ rho1 <- rho.IAPWS95(T=rep(this.T1, 2), P=c(p1, p2))
+ rho2 <- rho.IAPWS95(T=rep(this.T2, 2), P=c(p1, p2))
e1 <- water.AW90(T=rep(this.T1,2),rho=rho1,P=convert(c(p1,p2),'MPa'))
e2 <- water.AW90(T=rep(this.T2,2),rho=rho2,P=convert(c(p1,p2),'MPa'))
#p1 <- convert(-(1/e1[2]-1/e1[1])/(2*dp),'cm3bar')
@@ -379,22 +357,27 @@
### main loop; init dataframe output and density holders
w.out <- NULL
my.rho <- NULL
- # get densities and tell about it
- if(!quiet) msgout(" rho")
- my.rho <- rho()
+ # get densities unless only Psat is requested
+ if(!identical(tolower(property), "psat")) {
+ # calculate values of P for Psat
+ if(identical(P, "Psat")) P <- psat()
+ msgout(" rho")
+ my.rho <- rho.IAPWS95(T, P)
+ rho <- function() my.rho
+ }
for(i in 1:length(property)) {
- if(property[i] %in% c('e','kt')) {
- # expansivity isn't in the table yet... set it to zero
- warning('water: values of ',property[i],' are NA\n',call.=FALSE)
- inew <- rep(NA,length(T))
+ if(tolower(property[i]) %in% c("e", "kt")) {
+ # E and kT aren't here yet... set them to NA
+ warning("water.IAPWS95: values of ", property[i], " are NA\n", call.=FALSE)
+ inew <- rep(NA, length(T))
} else {
- if(!quiet) msgout(paste(" ", property[i], sep=""))
- inew <- get(property[i])()
+ msgout(paste(" ", property[i], sep=""))
+ inew <- get(tolower(property[i]))()
}
- #if(NA %in% inew) na.h2o <- TRUE
wnew <- data.frame(inew)
- if(i > 1) w.out <- cbind(w.out,wnew) else w.out <- wnew
+ if(i > 1) w.out <- cbind(w.out, wnew) else w.out <- wnew
}
- if(!quiet) msgout("\n")
+ msgout("\n")
+ names(w.out) <- property
return(w.out)
}
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/data/opt.csv 2013-02-12 22:35:25 UTC (rev 42)
@@ -1,2 +1,2 @@
Tr,Pr,Theta,Psi,R,cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol
-298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT,100,1,1
+298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT92,100,1,1
Modified: pkg/CHNOSZ/demo/CO2Ac.R
===================================================================
--- pkg/CHNOSZ/demo/CO2Ac.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/demo/CO2Ac.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -5,7 +5,7 @@
species("acetic acid", -3)
a <- affinity(O2=c(-85, -70, 4), T=c(25, 100, 4))
# hacking to write a title with formulas and subscripts
-lCO2 <- axis.label("CO2")[[1]]
+lCO2 <- axis.label("CO2")
main <- substitute(a~~b~~c,list(a=lCO2, b="buffered by",
c="acetic acid"))
d <- diagram(a, what="CO2", main=main)
@@ -14,7 +14,7 @@
d <- diagram(a, what="CO2", add=TRUE, lty=2)
# add a legend
lAC <- expr.species("CH3COOH", log="aq")
-ltext <- c(lAC, -3, -10)
+ltext <- c(as.expression(lAC), -3, -10)
lty <- c(NA, 1, 2)
legend("topright", legend=ltext, lty=lty, bg="white")
# do return.buffer and diagram(what) give the same results?
Modified: pkg/CHNOSZ/demo/nonideal.R
===================================================================
--- pkg/CHNOSZ/demo/nonideal.R 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/demo/nonideal.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -37,4 +37,4 @@
"charged species at 0, 25, 40 deg C, after Alberty, 2003",
sep="\n"),cex.main=0.95)
legend("topright", lty=c(NA, 1, 1, 1), col=c(NA, "blue", "black", "red"),
- legend=c(axis.label("T"), 0, 25, 40))
+ legend=c(as.expression(axis.label("T")), 0, 25, 40))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2013-02-06 15:17:21 UTC (rev 41)
+++ pkg/CHNOSZ/inst/NEWS 2013-02-12 22:35:25 UTC (rev 42)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9-9.4 (2013-02-06)
+CHANGES IN CHNOSZ 0.9-9.5 (2013-02-13)
--------------------------------------
- Fix calculation of free energy derivative in wjd().
@@ -36,6 +36,18 @@
- Examples of calculation of affinity of formation of CSG_METVO
(following Dick and Shock, 2011) added to protein.info.Rd.
+- Use consistent names for water properties (Speed, diel, QBorn, ...).
+
+- Add water.props() to get names of properties of water.
+
+- Remove 'isat' argument from water.SUPCRT92(); function now accepts
+ 'Psat' as value for 'P' argument.
+
+- Separate rho.IAPWS95() from water.IAPWS95().
+
+- Any property of water can be used as a variable in EOSvar().
+
+
CHANGES IN CHNOSZ 0.9-9 (2013-01-01)
------------------------------------
Added: pkg/CHNOSZ/inst/tests/test-EOSregress.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-EOSregress.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-EOSregress.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -0,0 +1,6 @@
+context("EOSregress")
+
+test_that("EOSvar stops with unknown variables", {
+ expect_error(EOSvar("TX", T=25, P=1), "can't find a variable named TX")
+})
+
Added: pkg/CHNOSZ/inst/tests/test-water.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-water.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-water.R 2013-02-12 22:35:25 UTC (rev 42)
@@ -0,0 +1,23 @@
+context("water")
+
+test_that("water.SUPCRT92() gives expected erros and warnings", {
+ expect_error(water.SUPCRT92("X"), "not available: X")
+ expect_error(water.SUPCRT92("Psat"), "please set P='Psat'")
+})
+
+test_that("water.SUPCRT92() gives expected values for E and kT", {
+ # E = V * alpha, kT = V * beta
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 42
More information about the CHNOSZ-commits
mailing list