[CHNOSZ-commits] r228 - in pkg/CHNOSZ: . R demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 29 12:58:24 CEST 2017
Author: jedick
Date: 2017-09-29 12:58:23 +0200 (Fri, 29 Sep 2017)
New Revision: 228
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/demo/DEW.R
Log:
code cleanup: remove unexported eos.args()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-29 10:58:23 UTC (rev 228)
@@ -1,6 +1,6 @@
Date: 2017-09-29
Package: CHNOSZ
-Version: 1.1.0-26
+Version: 1.1.0-27
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/R/cgl.R 2017-09-29 10:58:23 UTC (rev 228)
@@ -8,9 +8,6 @@
thermo <- get("thermo")
Tr <- thermo$opt$Tr
Pr <- thermo$opt$Pr
- eargs <- eos.args("mk", property = property)
- prop <- eargs$prop
- EOS.Prop <- eargs$Prop
# the number of T, P conditions
ncond <- max(c(length(T), length(P)))
@@ -22,13 +19,13 @@
# the parameters for *this* species
PAR <- parameters[k, ]
# start with NA values
- values <- data.frame(matrix(NA, ncol = length(prop), nrow=ncond))
- colnames(values) <- EOS.Prop
+ values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
+ colnames(values) <- property
# additional calculations for quartz and coesite
qtz <- quartz_coesite(PAR, T, P)
isqtz <- !identical(qtz$V, 0)
- for(i in 1:length(prop)) {
- PROP <- prop[i]
+ for(i in 1:length(property)) {
+ PROP <- property[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
@@ -36,20 +33,20 @@
# if at least one of the EoS parameters is available, zero out any NA's in the rest
if(hasEOS) PAR[, 13:19][, is.na(PAR[, 13:19])] <- 0
# equations for lambda adapted from HOK+98
- if(PROP == "cp") {
+ if(PROP == "Cp") {
# use constant Cp if the EoS parameters are not available
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
}
- if(PROP == "v") {
+ if(PROP == "V") {
if(isqtz) p <- qtz$V
else p <- rep(PAR$V, ncond)
}
- if(PROP %in% c("e", "kt")) {
+ if(PROP %in% c("E", "kT")) {
p <- rep(NA, ncond)
warning("cgl: E and/or kT of cr, gas and/or liq species are NA.")
}
- if(PROP == "g") {
+ if(PROP == "G") {
# use constant Cp if the EoS parameters are not available
if(!hasEOS) p <- PAR$Cp * (T - Tr - T * log(T/Tr)) else {
# Gibbs energy integral: the value at Tref plus heat capacity terms
@@ -70,7 +67,7 @@
else if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
p <- PAR$G + p
}
- if(PROP == "h") {
+ if(PROP == "H") {
# use constant Cp if the EoS parameters are not available
if(!hasEOS) p <- PAR$Cp * (T - Tr) else {
p <- PAR$a * (T - Tr) + PAR$b * (T^2 - Tr^2) / 2 +
@@ -86,7 +83,7 @@
#else p <- p + convert(PAR$V*(P-Pr),'calories')
p <- PAR$H + p
}
- if(PROP=="s") {
+ if(PROP=="S") {
# use constant Cp if the EoS parameters are not available
if(!hasEOS) p <- PAR$Cp * log(T/Tr) else {
p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-29 10:58:23 UTC (rev 228)
@@ -6,7 +6,7 @@
#source("util.args.R")
hkf <- function(property = NULL, parameters = NULL, T = 298.15, P = 1,
- contrib = c('n', 's', 'o'), H2O.props="rho") {
+ contrib = c("n", "s", "o"), H2O.props="rho") {
# calculate G, H, S, Cp, V, kT, and/or E using
# the revised HKF equations of state
# H2O.props - H2O properties needed for subcrt() output
@@ -16,19 +16,14 @@
Pr <- thermo$opt$Pr
Theta <- thermo$opt$Theta
Psi <- thermo$opt$Psi
- # argument handling
- eargs <- eos.args('hkf', property)
- property <- eargs$prop
- EOS.props <- eargs$props
- EOS.Props <- eargs$Prop
# make T and P equal length
if(!identical(P, "Psat")) {
if(length(P) < length(T)) P <- rep(P, length.out = length(T))
if(length(T) < length(P)) T <- rep(T, length.out = length(P))
}
# nonsolvation, solvation, and origination contribution
- notcontrib <- ! contrib %in% c('n', 's', 'o')
- if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o'); got", c2s(contrib[notcontrib])))
+ notcontrib <- ! contrib %in% c("n", "s", "o")
+ if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o); got", c2s(contrib[notcontrib])))
# get water properties
# rho - for subcrt() output and g function
# Born functions and epsilon - for HKF calculations
@@ -101,18 +96,18 @@
# various contributions to the properties
if(icontrib == "n") {
# nonsolvation ghs equations
- if(PROP == "h") {
+ 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") {
+ } 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)*(PAR$a3*(P-Pr)+PAR$a4*log((Psi+P)/(Psi+Pr)))
p <- p.c + p.a
- } else if(PROP == "g") {
+ } 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))) )
@@ -122,53 +117,53 @@
# at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA
if(!is.na(PAR$G)) p[T==Tr & P==Pr] <- 0
# nonsolvation cp v kt e equations
- } else if(PROP == 'cp') {
+ } 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')
+ } 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(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') *
- H2O.PT$QBorn + convert(dwdP,'cm3bar') * (-ZBorn - 1)
+ 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$NBorn
- if(PROP == 'e') p <- -convert(omega,'cm3bar') * H2O.PT$UBorn
+ # (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$NBorn
+ if(PROP == "E") p <- -convert(omega, "cm3bar") * H2O.PT$UBorn
}
- if(icontrib == 'o') {
+ if(icontrib == "o") {
# origination ghs equations
- if(PROP == 'g') {
+ if(PROP == "G") {
p <- PAR$G - PAR$S * (T-Tr)
- # don't inherit NA from PAR$S at Tr
+ # don"t inherit NA from PAR$S at Tr
p[T==Tr] <- PAR$G
}
- else if(PROP == 'h') p <- PAR$H
- else if(PROP == 's') p <- PAR$S
- # origination eos equations: senseless
- else if(PROP %in% tolower(EOS.props)) p <- 0 * T
+ else if(PROP == "H") p <- PAR$H
+ else if(PROP == "S") p <- PAR$S
+ # origination eos equations (Cp, V, kT, E): senseless
+ else p <- 0 * T
}
# accumulate the contribution
hkf.p <- hkf.p + p
@@ -176,7 +171,7 @@
wnew <- data.frame(hkf.p)
if(i > 1) w <- cbind(w, wnew) else w <- wnew
}
- colnames(w) <- EOS.Props
+ colnames(w) <- property
aq.out[[k]] <- w
}
return(list(aq=aq.out, H2O=H2O.PT))
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/R/subcrt.R 2017-09-29 10:58:23 UTC (rev 228)
@@ -9,9 +9,9 @@
#source("util.units.R")
#source("util.data.R")
-subcrt <- function(species, coeff=1, state=NULL, property=c('logK','G','H','S','V','Cp'),
- T=seq(273.15,623.15,25), P='Psat', grid=NULL, convert=TRUE, exceed.Ttr=FALSE,
- logact=NULL, action.unbalanced='warn', IS=0) {
+subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
+ T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
+ logact = NULL, action.unbalanced = "warn", IS = 0) {
# revise the call if the states have
# come as the second argument
@@ -47,12 +47,12 @@
}
# allowed properties
- properties <- c('rho','logK','G','H','S','Cp','V','kT','E')
+ properties <- c("rho", "logK", "G", "H", "S", "Cp", "V", "kT", "E")
# property checking
- outprop <- tolower(property)
- notproperty <- property[!outprop %in% tolower(properties)]
- if(length(notproperty) > 0) stop(paste(notproperty,
- 'are not valid properties\ntry rho, logK, G, H, S, V, Cp, kT, or E (or their lowercase equivalents)'))
+ calcprop <- property
+ notprop <- property[!calcprop %in% properties]
+ if(length(notprop) > 0) stop(paste(notprop,
+ "are not valid properties\ntry rho, logK, G, H, S, V, Cp, kT, or E"))
# length checking
if(do.reaction & length(species)!=length(coeff))
stop('coeff must be same length as the number of species.')
@@ -262,14 +262,14 @@
# calculate the properties
# if we want affinities we must have logK; include it in the ouput
- if(!is.null(logact)) if(!'logk' %in% outprop) outprop <- c('logk', outprop)
- # if logK but not g was requested, we need to calculate g
- eosprop <- outprop
- if('logk' %in% outprop & ! 'g' %in% outprop) eosprop <- c(eosprop, 'g')
+ if(!is.null(logact)) if(!'logK' %in% calcprop) calcprop <- c('logK', calcprop)
+ # if logK but not G was requested, we need to calculate G
+ eosprop <- calcprop
+ if('logK' %in% calcprop & ! 'G' %in% calcprop) eosprop <- c(eosprop, 'G')
# also get g if we are dealing with mineral phases
- if(!'g' %in% eosprop & length(inpho) > length(sinfo)) eosprop <- c(eosprop, 'g')
- # don't request logk or rho from the eos ...
- eosprop <- eosprop[!eosprop %in% c('logk','rho')]
+ if(!'G' %in% eosprop & length(inpho) > length(sinfo)) eosprop <- c(eosprop, 'G')
+ # don't request logK or rho from the eos ...
+ eosprop <- eosprop[!eosprop %in% c('logK','rho')]
# the reaction result will go here
out <- list()
# aqueous species and H2O properties
@@ -302,7 +302,7 @@
p.cgl <- cgl(eosprop, parameters = si, T = T, P = P)
# replace Gibbs energies with NA where the
# phases are beyond their temperature range
- if('g' %in% eosprop) {
+ if('G' %in% eosprop) {
# 20080304 this code is weird and hard to read - needs a lot of cleanup!
# 20120219 cleaned up somewhat; using exceed.Ttr and NA instead of do.phases and 999999
# the numbers of the cgl species (becomes 0 for any that aren't cgl)
@@ -350,7 +350,7 @@
# water
if(any(isH2O)) {
- p.H2O <- H2O.PT[, match(eosprop, tolower(colnames(H2O.PT))), drop=FALSE]
+ p.H2O <- H2O.PT[, match(eosprop, colnames(H2O.PT)), drop=FALSE]
p.H2O <- list(p.H2O)
out <- c(out, rep(p.H2O, sum(isH2O == TRUE)))
}
@@ -362,7 +362,7 @@
}
# logK
- if('logk' %in% outprop) {
+ if('logK' %in% calcprop) {
for(i in 1:length(out)) {
# NOTE: the following depends on the water function renaming g to G
out[[i]] <- cbind(out[[i]],data.frame(logK=convert(out[[i]]$G,'logK',T=T)))
@@ -459,9 +459,9 @@
isH2O <- isH2O.new
# the order of the properties
- if(length(outprop)>1) for(i in 1:length(out)) {
+ if(length(calcprop)>1) for(i in 1:length(out)) {
# keep state/loggam columns if they exists
- ipp <- match(outprop,tolower(colnames(out[[i]])))
+ ipp <- match(calcprop, colnames(out[[i]]))
if('state' %in% colnames(out[[i]])) ipp <- c(ipp,match('state',colnames(out[[i]])))
if('loggam' %in% colnames(out[[i]])) ipp <- c(ipp,match('loggam',colnames(out[[i]])))
out[[i]] <- out[[i]][,ipp,drop=FALSE]
@@ -524,7 +524,7 @@
if(convert) P.out <- outvert(P,"bar") else P.out <- P
# try to stuff in a column of rho if we have aqueous species
# watch out! supcrt-ish densities are in g/cc not kg/m3
- if('rho' %in% outprop | (missing(property) & any(c(isaq,isH2O))) & (names(out)[i])!='state')
+ if('rho' %in% calcprop | (missing(property) & any(c(isaq,isH2O))) & (names(out)[i])!='state')
out[[i]] <- cbind(data.frame(T=T.out,P=P.out,rho=H2O.PT$rho/1000),out[[i]])
else
out[[i]] <- cbind(data.frame(T=T.out,P=P.out,out[[i]]))
Modified: pkg/CHNOSZ/R/util.args.R
===================================================================
--- pkg/CHNOSZ/R/util.args.R 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/R/util.args.R 2017-09-29 10:58:23 UTC (rev 228)
@@ -3,44 +3,6 @@
### unexported functions ###
-# These functions are used to normalize user-input arguments, which are case-insensitive.
-
-# return a list with elements:
-# `prop` for all the properties available for the specified equations-of-state
-# `prop` for the lower-case version of property
-# `Prop` for the upper-case (of first letter) version of property
-# produces an error if any of `property` is not in the list of available properties.
-# (See water() and subcrt() for the available properties for different species.)
-eos.args <- function(eos='',property=NULL,T=NULL,P=NULL) {
- # the available properties for supcrt, probably
- props <- c('G','H','S','Cp','V','kT','E')
- if(eos=='water') {
- # things we also get with water
- props <- c(props,'A','U','Cv','Psat','rho','Q','X','Y','epsilon','w')
- # they keep on coming: things we also get with SUPCRT92
- if(get("thermo")$opt$water == "SUPCRT92")
- props <- c(props,'Z','visc','tcond','tdiff','Prndtl','visck','albe','daldT','alpha','beta')
- else
- props <- c(props,'P','N','UBorn','de.dT','de.dP')
- }
- # so others can find out what we're about
- if(is.null(property)) return(props)
- # default: all properties and contributions
- if(is.null(property)) property <- props
- # the lowercase equivalent of the argument
- prop <- tolower(property)
- # and returns its lower- and upper- case equivalents
- # (prop and Prop) and the available properties
- Prop <- props[match(prop,tolower(props))]
- #contrib <- tolower(contrib)
- # error: not a property
- notprop <- ! prop %in% tolower(props)
- if(TRUE %in% notprop)
- stop('thermo.args: properties ',c2s(prop[notprop]),' not in ',c2s(props),'\n')
- # return arguments
- return(list(props=props,prop=prop,Prop=Prop))
-}
-
# force T and P to equal length
# also looks for the keyword Psat in the value of P and substitutes calculated values of the saturation vapor pressure
TP.args <- function(T=NULL, P=NULL) {
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2017-09-29 09:47:24 UTC (rev 227)
+++ pkg/CHNOSZ/demo/DEW.R 2017-09-29 10:58:23 UTC (rev 228)
@@ -41,7 +41,7 @@
text(lastT+25, tail(logm, 1), Pkb, adj=0)
}
t1 <- quote("Solubility of"~alpha*"-quartz")
-t2 <- "(Sverjensky et al., 2014)"
+t2 <- "(after Sverjensky et al., 2014)"
mtitle(as.expression(c(t1, t2)))
# TODO: lines are a little low at highest P and P ...
# does the Berman, 1988 quartz data increase high-PT solubilities?
More information about the CHNOSZ-commits
mailing list