[CHNOSZ-commits] r210 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 25 05:44:27 CEST 2017
Author: jedick
Date: 2017-09-25 05:44:26 +0200 (Mon, 25 Sep 2017)
New Revision: 210
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/tests/testthat/test-species.R
Log:
hkf() and cgl(): 'ghs' and 'eos' arguments merged into 'parameters'
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-25 03:44:26 UTC (rev 210)
@@ -1,6 +1,6 @@
Date: 2017-09-24
Package: CHNOSZ
-Version: 1.1.0-8
+Version: 1.1.0-9
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/R/EOSregress.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -5,14 +5,14 @@
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) (calories)
- Cp_s <- hkf("Cp", T=T, P=P, eos=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ Cp_s <- hkf("Cp", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
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", T=T, P=P, eos=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ V_s <- hkf("V", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
return(V_s[[1]][, 1]/convert(omega.PrTr, "cm3bar"))
}
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/R/cgl.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -2,7 +2,7 @@
# calculate standard thermodynamic properties of non-aqueous species
# 20060729 jmd
-cgl <- function(property=NULL,T=298.15,P=1,ghs=NULL,eos=NULL) {
+cgl <- function(property=NULL,T=298.15,P=1,parameters=NULL) {
# calculate properties of crystalline, liquid (except H2O) and gas species
# argument handling
thermo <- get("thermo")
@@ -10,79 +10,77 @@
Pr <- thermo$opt$Pr
eargs <- eos.args('mk',property=property)
prop <- eargs$prop
- props <- eargs$props
- Prop <- eargs$Prop
+ EOS.Prop <- eargs$Prop
# a list - the result
x <- list()
- for(k in 1:nrow(ghs)) {
+ for(k in 1:nrow(parameters)) {
# loop over each species
- GHS <- ghs[k,]
- EOS <- eos[k,]
+ PAR <- parameters[k,]
w <- NULL
for(i in 1:length(prop)) {
- property <- prop[i]
+ PROP <- prop[i]
# a test for availability of the EoS parameters
# here we assume that the parameters are in the same position as in thermo$obigt
# leave T transition (in 20th column) alone
- hasEOS <- any(!is.na(EOS[, 13:19]))
+ hasEOS <- any(!is.na(PAR[, 13:19]))
# if at least one of the EoS parameters is available, zero out any NA's in the rest
- if(hasEOS) EOS[, 13:19][, is.na(EOS[, 13:19])] <- 0
+ if(hasEOS) PAR[, 13:19][, is.na(PAR[, 13:19])] <- 0
# equations for lambda adapted from HOK+98
- if(property=='cp') {
+ if(PROP=='cp') {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- EOS$Cp
- else p <- EOS$a + EOS$b*T + EOS$c*T^-2 + EOS$d*T^-0.5 + EOS$e*T^2 + EOS$f*T^EOS$lambda
- } else if(property=='v') {
- p <- rep(EOS$V,length(T))
- } else if(property %in% c('e','kt')) {
+ if(!hasEOS) p <- PAR$Cp
+ else p <- PAR$a + PAR$b*T + PAR$c*T^-2 + PAR$d*T^-0.5 + PAR$e*T^2 + PAR$f*T^PAR$lambda
+ } else if(PROP=='v') {
+ p <- rep(PAR$V,length(T))
+ } else if(PROP %in% c('e','kt')) {
p <- rep(NA,length(T))
warning('cgl: E and/or kT of cr, gas and/or liq species are NA.')
- } else if(property=='g') {
+ } else if(PROP=='g') {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- GHS$G + EOS$Cp*(T-Tr-T*log(T/Tr)) else
+ if(!hasEOS) p <- PAR$G + PAR$Cp*(T-Tr-T*log(T/Tr)) else
# Gibbs energy integral: the value at Tref plus heat capacity terms
- p <- GHS$G + EOS$a*(T-Tr-T*log(T/Tr)) -
- EOS$b*(T-Tr)^2/2 - EOS$c*(1/T + T/Tr^2 - 2/Tr)/2 -
- EOS$d*(T^0.5-0.5*T*Tr^-0.5-0.5*Tr^0.5)/-0.25 -
- EOS$e*(T^3-3*T*Tr^2+2*Tr^3)/6
+ p <- PAR$G + PAR$a*(T-Tr-T*log(T/Tr)) -
+ PAR$b*(T-Tr)^2/2 - PAR$c*(1/T + T/Tr^2 - 2/Tr)/2 -
+ PAR$d*(T^0.5-0.5*T*Tr^-0.5-0.5*Tr^0.5)/-0.25 -
+ PAR$e*(T^3-3*T*Tr^2+2*Tr^3)/6
# entropy and volume terms
- if(!is.na(GHS$S)) p <- p - GHS$S*(T-Tr)
- if(!is.na(EOS$V)) p <- p + convert(EOS$V*(P-Pr),'calories')
+ if(!is.na(PAR$S)) p <- p - PAR$S*(T-Tr)
+ if(!is.na(PAR$V)) p <- p + convert(PAR$V*(P-Pr),'calories')
# use additional heat capacity term if it's defined
- if(!is.na(EOS$f) & !is.na(EOS$lambda)) if(EOS$f!=0) {
- if(EOS$lambda== -1) p <- p + EOS$f*(log(T/Tr)-T*(1/Tr-1/T))
- else p <- p + EOS$f * ( T^(EOS$lambda+1) - (EOS$lambda+1)*T*Tr^EOS$lambda +
- EOS$lambda*Tr^(EOS$lambda+1) ) / ( EOS$lambda*(EOS$lambda+1) )
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f!=0) {
+ if(PAR$lambda== -1) p <- p + PAR$f*(log(T/Tr)-T*(1/Tr-1/T))
+ else p <- p + PAR$f * ( T^(PAR$lambda+1) - (PAR$lambda+1)*T*Tr^PAR$lambda +
+ PAR$lambda*Tr^(PAR$lambda+1) ) / ( PAR$lambda*(PAR$lambda+1) )
}
- } else if(property=='h') {
+ } else if(PROP=='h') {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- EOS$Cp*(T-Tr) else
- p <- EOS$a*(T-Tr) + EOS$b*(T^2-Tr^2)/2 +
- EOS$c*(1/T-1/Tr)/-1 + EOS$d*(T^0.5-Tr^0.5)/0.5 +
- EOS$e*(T^3-Tr^3)/3
+ if(!hasEOS) p <- PAR$Cp*(T-Tr) else
+ p <- PAR$a*(T-Tr) + PAR$b*(T^2-Tr^2)/2 +
+ PAR$c*(1/T-1/Tr)/-1 + PAR$d*(T^0.5-Tr^0.5)/0.5 +
+ PAR$e*(T^3-Tr^3)/3
# SUPCRT seems to ignore this term? ... 20070802
- # + convert(EOS$V*(P-Pr),'calories')
- p <- GHS$H + p
- if(!is.na(EOS$f) & !is.na(EOS$lambda)) if(EOS$f!=0) {
- if(EOS$lambda == -1) p <- p + EOS$f*log(T/Tr)
- else p <- p - EOS$f * ( T^(EOS$lambda+1) - Tr^(EOS$lambda+1) ) / (EOS$lambda+1)
+ # + convert(PAR$V*(P-Pr),'calories')
+ p <- PAR$H + p
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f!=0) {
+ if(PAR$lambda == -1) p <- p + PAR$f*log(T/Tr)
+ else p <- p - PAR$f * ( T^(PAR$lambda+1) - Tr^(PAR$lambda+1) ) / (PAR$lambda+1)
}
- } else if(property=='s') {
+ } else if(PROP=='s') {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- EOS$Cp*log(T/Tr) else
- p <- EOS$a*log(T/Tr) + EOS$b*(T-Tr) +
- EOS$c*(T^(-2)-Tr^(-2))/(-2) + EOS$e*(T^2-Tr^2)/2 +
- EOS$d*(T^-0.5-Tr^-0.5)/-0.5
- p <- GHS$S + p
- if(!is.na(EOS$f) & !is.na(EOS$lambda)) if(EOS$f!=0) {
- p <- p + EOS$f*(T^EOS$lambda-Tr^EOS$lambda)/EOS$lambda
+ if(!hasEOS) p <- PAR$Cp*log(T/Tr) else
+ p <- PAR$a*log(T/Tr) + PAR$b*(T-Tr) +
+ PAR$c*(T^(-2)-Tr^(-2))/(-2) + PAR$e*(T^2-Tr^2)/2 +
+ PAR$d*(T^-0.5-Tr^-0.5)/-0.5
+ p <- PAR$S + p
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f!=0) {
+ p <- p + PAR$f*(T^PAR$lambda-Tr^PAR$lambda)/PAR$lambda
}
}
wnew <- data.frame(p)
if(i>1) w <- cbind(w,wnew) else w <- wnew
}
- colnames(w) <- Prop
+ colnames(w) <- EOS.Prop
x[[k]] <- w
}
return(x)
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -2,7 +2,7 @@
# calculate thermodynamic properties using equations of state
# 11/17/03 jmd
-hkf <- function(property = NULL, T = 298.15, P = 1, ghs = NULL, eos = NULL,
+hkf <- function(property = NULL, T = 298.15, P = 1, parameters = NULL,
contrib = c('n', 's', 'o'), H2O.PT = NULL, H2O.PrTr = NULL, domega = TRUE) {
# calculate G, H, S, Cp, V, kT, and/or E using
# the revised HKF equations of state
@@ -15,25 +15,25 @@
# argument handling
eargs <- eos.args('hkf', property)
property <- eargs$prop
- props <- eargs$props
- Prop <- eargs$Prop
- domega <- rep(domega, length.out = nrow(eos))
+ EOS.props <- eargs$props
+ EOS.Props <- eargs$Prop
+ domega <- rep(domega, length.out = nrow(parameters))
# nonsolvation, solvation, and origination contribution
- contribs <- c('n', 's', 'o')
- notcontrib <- ! contrib %in% contribs
+ notcontrib <- ! contrib %in% c('n', 's', 'o')
if(TRUE %in% notcontrib)
- stop(paste('argument', c2s(contrib[notcontrib]), 'not in', c2s(contribs), 'n'))
+ stop(paste("contrib must be in c('n', 's', 'o'); got", c2s(contrib[notcontrib])))
# get water properties, if they weren't supplied in arguments (and we want solvation props)
+ dosupcrt <- thermo$opt$water != "IAPWS95"
if('s' %in% contrib) {
H2O.props <- c("QBorn", "XBorn", "YBorn", "diel")
- # only take these ones if we're in SUPCRT92 compatibility mode
- dosupcrt <- thermo$opt$water != "IAPWS95"
if(dosupcrt) {
+ # using H2O92D.f from SUPCRT92
# (rho, alpha, daldT, beta - for partial derivatives of omega (g function))
H2O.props <- c(H2O.props, "rho", "alpha", "daldT", "beta")
} else {
+ # using IAPWS-95
# (NBorn, UBorn - for compressibility, expansibility)
- H2O.props <- c(H2O.props,'NBorn','UBorn')
+ 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)
@@ -42,34 +42,31 @@
}
# a list to store the result
x <- list()
- nspecies <- nrow(ghs)
- # we can be called with NULL ghs (by Cp_s_var, V_s_var)
- if(is.null(nspecies)) nspecies <- nrow(eos)
+ nspecies <- nrow(parameters)
for(k in 1:nspecies) {
# loop over each species
- GHS <- ghs[k, ]
- EOS <- eos[k, ]
+ PAR <- parameters[k, ]
# substitute Cp and V for missing EoS parameters
# here we assume that the parameters are in the same position as in thermo$obigt
# 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
+ if(all(is.na(PAR[, 17:18]))) PAR[, 17] <- PAR$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")
+ if(all(is.na(PAR[, 13:16]))) PAR[, 13] <- convert(PAR$V, "calories")
# test for availability of the EoS parameters
- hasEOS <- any(!is.na(EOS[, 13:20]))
+ hasEOS <- any(!is.na(PAR[, 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
+ if(hasEOS) PAR[, 13:20][, is.na(PAR[, 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
+ omega <- PAR$omega # omega.PrTr
# its derivatives are zero unless the g function kicks in
dwdP <- dwdT <- d2wdT2 <- numeric(length(T))
- Z <- EOS$Z
- omega.PT <- rep(EOS$omega, length(T))
- if(!is.na(Z)) if(Z != 0) if(domega[k]) if(dosupcrt) {
+ Z <- PAR$Z
+ omega.PT <- rep(PAR$omega, length(T))
+ if(!identical(Z, 0) & domega[k] & dosupcrt) {
# g and f function stuff (Shock et al., 1992; Johnson et al., 1992)
rhohat <- H2O.PT$rho/1000 # just converting kg/m3 to g/cm3
g <- gfun(rhohat, convert(T, "C"), P, H2O.PT$alpha, H2O.PT$daldT, H2O.PT$beta)
@@ -87,81 +84,81 @@
# loop over each property
w <- NULL
for(i in 1:length(property)) {
- prop <- property[i]
+ PROP <- property[i]
# over nonsolvation, solvation, or origination contributions
hkf.p <- numeric(length(T))
for(icontrib in contrib) {
# various contributions to the properties
if(icontrib == "n") {
# nonsolvation ghs equations
- if(prop == "h") {
- p.c <- EOS$c1*(T-Tr) - EOS$c2*(1/(T-Theta)-1/(Tr-Theta))
- p.a <- EOS$a1*(P-Pr) + EOS$a2*log((Psi+P)/(Psi+Pr)) +
- ((2*T-Theta)/(T-Theta)^2)*(EOS$a3*(P-Pr)+EOS$a4*log((Psi+P)/(Psi+Pr)))
+ if(PROP == "h") {
+ p.c <- PAR$c1*(T-Tr) - PAR$c2*(1/(T-Theta)-1/(Tr-Theta))
+ p.a <- PAR$a1*(P-Pr) + PAR$a2*log((Psi+P)/(Psi+Pr)) +
+ ((2*T-Theta)/(T-Theta)^2)*(PAR$a3*(P-Pr)+PAR$a4*log((Psi+P)/(Psi+Pr)))
p <- p.c + p.a
- } else if(prop == "s") {
- p.c <- EOS$c1*log(T/Tr) -
- (EOS$c2/Theta)*( 1/(T-Theta)-1/(Tr-Theta) +
+ } else if(PROP == "s") {
+ p.c <- PAR$c1*log(T/Tr) -
+ (PAR$c2/Theta)*( 1/(T-Theta)-1/(Tr-Theta) +
log( (Tr*(T-Theta))/(T*(Tr-Theta)) )/Theta )
- p.a <- (T-Theta)^(-2)*(EOS$a3*(P-Pr)+EOS$a4*log((Psi+P)/(Psi+Pr)))
+ p.a <- (T-Theta)^(-2)*(PAR$a3*(P-Pr)+PAR$a4*log((Psi+P)/(Psi+Pr)))
p <- p.c + p.a
- } else if(prop == "g") {
- p.c <- -EOS$c1*(T*log(T/Tr)-T+Tr) -
- EOS$c2*( (1/(T-Theta)-1/(Tr-Theta))*((Theta-T)/Theta) -
+ } else if(PROP == "g") {
+ p.c <- -PAR$c1*(T*log(T/Tr)-T+Tr) -
+ PAR$c2*( (1/(T-Theta)-1/(Tr-Theta))*((Theta-T)/Theta) -
(T/Theta^2)*log((Tr*(T-Theta))/(T*(Tr-Theta))) )
- p.a <- EOS$a1*(P-Pr) + EOS$a2*log((Psi+P)/(Psi+Pr)) +
- (EOS$a3*(P-Pr) + EOS$a4*log((Psi+P)/(Psi+Pr)))/(T-Theta)
+ p.a <- PAR$a1*(P-Pr) + PAR$a2*log((Psi+P)/(Psi+Pr)) +
+ (PAR$a3*(P-Pr) + PAR$a4*log((Psi+P)/(Psi+Pr)))/(T-Theta)
p <- p.c + p.a
# 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(!is.na(PAR$G)) p[T==Tr & P==Pr] <- 0
# nonsolvation cp v kt e equations
- } else if(prop == 'cp') {
- p <- EOS$c1 + EOS$c2 * ( T - Theta ) ^ (-2)
- } else if(prop == 'v') {
- p <- convert(EOS$a1,'cm3bar') +
- convert(EOS$a2,'cm3bar') / ( Psi + P) +
- ( convert(EOS$a3,'cm3bar') + convert(EOS$a4,'cm3bar') / ( Psi + P ) ) / ( T - Theta)
- } else if(prop == 'kt') {
- p <- ( convert(EOS$a2,'cm3bar') +
- convert(EOS$a4,'cm3bar') / (T - Theta) ) * (Psi + P) ^ (-2)
- } else if(prop == 'e') {
- p <- convert( - ( EOS$a3 + EOS$a4 / convert((Psi + P),'calories') ) *
- (T - Theta) ^ (-2),'cm3bar')
+ } else if(PROP == 'cp') {
+ p <- PAR$c1 + PAR$c2 * ( T - Theta ) ^ (-2)
+ } else if(PROP == 'v') {
+ p <- convert(PAR$a1, 'cm3bar') +
+ convert(PAR$a2, 'cm3bar') / (Psi + P) +
+ (convert(PAR$a3, 'cm3bar') + convert(PAR$a4,'cm3bar') / (Psi + P)) / (T - Theta)
+ } else if(PROP == 'kt') {
+ p <- (convert(PAR$a2, 'cm3bar') +
+ convert(PAR$a4, 'cm3bar') / (T - Theta)) * (Psi + P) ^ (-2)
+ } else if(PROP == 'e') {
+ p <- convert( - (PAR$a3 + PAR$a4 / convert((Psi + P),'calories')) *
+ (T - Theta) ^ (-2), 'cm3bar')
}
}
if(icontrib == "s") {
# solvation ghs equations
- if(prop == "g") {
+ if(PROP == "g") {
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(!is.na(PAR$G)) p[T==Tr & P==Pr] <- 0
}
- if(prop == "h")
+ if(PROP == "h")
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")
+ if(PROP == "s")
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$XBorn + 2*T*H2O.PT$YBorn*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') *
+ if(PROP == 'v') p <- -convert(omega.PT,'cm3bar') *
H2O.PT$QBorn + convert(dwdP,'cm3bar') * (-ZBorn - 1)
# TODO: 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
- if(prop == 'e') p <- -convert(omega,'cm3bar') * H2O.PT$UBorn
+ if(PROP == 'kt') p <- convert(omega,'cm3bar') * H2O.PT$N
+ if(PROP == 'e') p <- -convert(omega,'cm3bar') * H2O.PT$UBorn
}
if(icontrib == 'o') {
# origination ghs equations
- if(prop == 'g') {
- p <- GHS$G - GHS$S * (T-Tr)
- # don't inherit NA from GHS$S at Tr
- p[T==Tr] <- GHS$G
+ if(PROP == 'g') {
+ p <- PAR$G - PAR$S * (T-Tr)
+ # don't inherit NA from PAR$S at Tr
+ p[T==Tr] <- PAR$G
}
- else if(prop == 'h') p <- GHS$H
- else if(prop == 's') p <- GHS$S
+ else if(PROP == 'h') p <- PAR$H
+ else if(PROP == 's') p <- PAR$S
# origination eos equations: senseless
- else if(prop %in% tolower(props)) p <- 0 * T
+ else if(PROP %in% tolower(EOS.props)) p <- 0 * T
}
# accumulate the contribution
hkf.p <- hkf.p + p
@@ -169,7 +166,7 @@
wnew <- data.frame(hkf.p)
if(i > 1) w <- cbind(w, wnew) else w <- wnew
}
- colnames(w) <- Prop
+ colnames(w) <- EOS.Props
x[[k]] <- w
}
return(x)
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/R/subcrt.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -288,7 +288,7 @@
#si <- info(inpho[isaq],quiet=TRUE)
si <- obigt2eos(thermo$obigt[inpho[isaq],], "aq", fixGHS = TRUE)
domega <- thermo$obigt$name[inpho[isaq]] != 'H+'
- p.aq <- hkf(eosprop,T=T,P=P,ghs=si,eos=si,H2O.PT=H2O.PT,H2O.PrTr=H2O.PrTr,domega=domega)
+ p.aq <- hkf(eosprop,T=T,P=P,parameters=si,H2O.PT=H2O.PT,H2O.PrTr=H2O.PrTr,domega=domega)
if(any(IS!=0)) p.aq <- nonideal(inpho[isaq],p.aq,newIS,T)
out <- c(out,p.aq)
}
@@ -300,7 +300,7 @@
if(TRUE %in% iscgl) {
#si <- info(inpho[iscgl],quiet=TRUE)
si <- obigt2eos(thermo$obigt[inpho[iscgl],], "cgl", fixGHS = TRUE)
- p.cgl <- cgl(eosprop,T=T,P=P,ghs=si,eos=si)
+ p.cgl <- cgl(eosprop, T = T, P = P, parameters = si)
# replace Gibbs energies with NA where the
# phases are beyond their temperature range
if('g' %in% eosprop) {
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/R/util.misc.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -48,18 +48,18 @@
# set the starting GHS to 0 (in case they're NA - we only need the increments over temperature)
thisinfo$G <- thisinfo$H <- thisinfo$S <- 0
# the HS increments from 298.15 to Ttr
- HSinc <- cgl(c("H", "S"), T=c(298.15, Ttr), ghs=thisinfo, eos=thisinfo)
+ HSinc <- cgl(c("H", "S"), T=c(298.15, Ttr), parameters=thisinfo)
Hf_Tr <- c(Hf_Tr, Hf - diff(HSinc[[1]]$H))
S_Tr <- c(S_Tr, S - diff(HSinc[[1]]$S))
# plug in the calculated S_Tr to calculate the G increment correctly
thisinfo$S <- tail(S_Tr, 1)
- Ginc <- cgl("G", T=c(298.15, Ttr), ghs=thisinfo, eos=thisinfo)
+ Ginc <- cgl("G", T=c(298.15, Ttr), parameters=thisinfo)
Gf_Tr <- c(Gf_Tr, Gf - diff(Ginc[[1]]$G))
}
# the temperature of the next transition
Ttr <- thisinfo$T
# the GHS increments from Tprev to Ttr
- GHCinc <- cgl(c("G", "H", "S"), T=c(Tprev, Ttr), ghs=thisinfo, eos=thisinfo)
+ GHCinc <- cgl(c("G", "H", "S"), T=c(Tprev, Ttr), parameters=thisinfo)
# the GHS + transition at Tr
Gf <- Gf + diff(GHCinc[[1]]$G)
Hf <- Hf + diff(GHCinc[[1]]$H) + Htr[i]
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-25 03:44:26 UTC (rev 210)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-8 (2017-09-24)
+CHANGES IN CHNOSZ 1.1.0-9 (2017-09-25)
--------------------------------------
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
@@ -21,6 +21,9 @@
inconsistent results for metastable steam (Zavarin et al., 2016,
LLNL-TR-701407, doi: 10.2172/1325873).
+- In hkf() and cgl(), combine 'ghs' and 'eos' arguments into single
+ argument named 'parameters'.
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/man/eos.Rd 2017-09-25 03:44:26 UTC (rev 210)
@@ -8,8 +8,8 @@
}
\usage{
- cgl(property = NULL, T = 298.15, P = 1, ghs = NULL, eos = NULL)
- hkf(property = NULL, T = 298.15, P = 1, ghs = NULL, eos = NULL,
+ cgl(property = NULL, T = 298.15, P = 1, parameters = NULL)
+ hkf(property = NULL, T = 298.15, P = 1, parameters = NULL,
contrib = c("n","s","o"), H2O.PT = NULL, H2O.PrTr = NULL,
domega = TRUE)
}
@@ -18,8 +18,7 @@
\item{property}{character, name(s) of properties to calculate}
\item{T}{numeric, temperature(s) at which to calculate properties (K)}
\item{P}{numeric, pressure(s) at which to calculate properties (bar)}
- \item{ghs}{dataframe, values of the standard molal Gibbs energy and enthalpy of formation from the elements and entropy at 25 \eqn{^{\circ}}{°}C and 1 bar}
- \item{eos}{dataframe, values of the equations-of-state parameters}
+ \item{parameters}{dataframe, standard molal thermodynamic properties and equations-of-state parameters in the format of \code{thermo$obigt}}
\item{contrib}{character, which contributions to consider in the revised HKF equations equations of state: (\code{n})onsolvation, (\code{s})olvation (the \eqn{\omega}{omega} terms), or (o)rigination contributions (i.e., the property itself at 25 \eqn{^{\circ}}{°}C and 1 bar). Default is \code{c("n","s","o")}, for all contributions}
\item{H2O.PT}{dataframe, values of the electrostatic properties of water at the temperature(s) and pressure(s) of interest}
\item{H2O.PrTr}{dataframe, values of the electrostatic properties of water at the reference temperature and pressure}
@@ -38,7 +37,7 @@
\code{H2O.PT} and \code{H2O.PrTr} are optional arguments that contain the electrostatic properties of \eqn{\mathrm{H_2O}}{H2O} required for the calculations.
If either of these is \code{NULL} (the default), the required values are retrieved using \code{\link{water}}.
-Unless \code{domega}, the value of which is recycled to the number of species (rows of ghs and eos), is FALSE for any of the species, the temperature and pressure derivatives of the \code{omega} parameter for charged species (where \code{Z != 0}) are calculated using the \eqn{g}{g}- and \eqn{f}{f}-functions described by Shock et al., 1992 and Johnson et al., 1992.
+Unless \code{domega}, the value of which is recycled to the number of species (number of rows of \code{parameters}), is FALSE for any of the species, the temperature and pressure derivatives of the \code{omega} parameter for charged species (where \code{Z != 0}) are calculated using the \eqn{g}{g}- and \eqn{f}{f}-functions described by Shock et al., 1992 and Johnson et al., 1992.
This option is turned off when the IAPWS-95 equations are activated (see \code{\link{water}}).
The parameters in the \code{cgl} equations of state for crystalline, gas and liquid species (except liquid water) include \code{V}, \code{a}, \code{b}, \code{c}, \code{d}, \code{e}, \code{f} and \code{lambda}.
@@ -59,7 +58,7 @@
}
\value{
-A list of length equal to the number of species (i.e., number rows of supplied \code{ghs} and \code{eos} values).
+A list of length equal to the number of species (i.e., number rows of \code{parameters}).
Each element of the list contains a dataframe, each column of which corresponds to one of the specified properties; the number of rows is equal to the number of pressure-temperature points.
}
@@ -70,22 +69,23 @@
\examples{
\dontshow{data(thermo)}
## aqueous species
-a <- info(info("methane","aq"))
-hkf(property="Cp",ghs=a,eos=a)
+CH4aq <- info(info("methane", "aq"))
+hkf(property = "Cp", parameters = CH4aq)
# the non-solvation heat capacity
-hkf(property="Cp",ghs=a,eos=a,contrib="n")
+hkf(property = "Cp", parameters = CH4aq, contrib = "n")
# at different temperature and pressure
-hkf(property="Cp",ghs=a,eos=a,T=c(373.15,473.15),P=1000)
+hkf(property = "Cp", parameters = CH4aq, T = c(373.15,473.15), P = 1000)
## crystalline, gas, liquid species
-a <- info(info("methane","gas"))
-cgl(property="Cp",ghs=a,eos=a)
+CH4gas <- info(info("methane", "gas"))
+cgl(property = "Cp", parameters = CH4gas)
# melting and vaporization of n-octane
-a <- info(info(rep("n-octane",3),c("cr","liq","gas")))
-b <- cgl(property="G",ghs=a,eos=a,T=seq(200,420,10),P=1)
-which.pmax(b,pmin=TRUE) # 1 = cr, 2 = liq, 3 = gas
+C8H18par <- info(info(rep("n-octane", 3), c("cr", "liq", "gas")))
+myT <- seq(200, 420, 10)
+DG0f <- cgl(property = "G", parameters = C8H18par, T = myT, P = 1)
+cbind(T = myT, which.pmax(DG0f, pmin = TRUE)) # 1 = cr, 2 = liq, 3 = gas
# compare that result with the tabulated transition temperatures
-print(a)
+print(C8H18par)
}
\references{
Modified: pkg/CHNOSZ/tests/testthat/test-species.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-species.R 2017-09-24 13:46:22 UTC (rev 209)
+++ pkg/CHNOSZ/tests/testthat/test-species.R 2017-09-25 03:44:26 UTC (rev 210)
@@ -60,5 +60,5 @@
# here it's "touched" (but not added or modified)
expect_equal(species("CO2", index.return=TRUE), 1)
expect_equal(species(c("H2O", "NH3"), index.return=TRUE), c(2, 3))
- expect_equal(species(2, "gas", index.return=TRUE), 2)
+ expect_equal(species(1, "gas", index.return=TRUE), 1)
})
More information about the CHNOSZ-commits
mailing list