[CHNOSZ-commits] r89 - in pkg/CHNOSZ: . R demo inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 13 19:27:57 CEST 2015
Author: jedick
Date: 2015-06-13 19:27:56 +0200 (Sat, 13 Jun 2015)
New Revision: 89
Added:
pkg/CHNOSZ/demo/ORP.R
pkg/CHNOSZ/inst/tests/test-mosaic.R
pkg/CHNOSZ/inst/tests/test-util.list.R
Removed:
pkg/CHNOSZ/demo/orp.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/iprotein.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/util.list.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/density.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-iprotein.R
pkg/CHNOSZ/man/examples.Rd
Log:
fix mosaic() to scale affinities by relative abundances of basis species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/DESCRIPTION 2015-06-13 17:27:56 UTC (rev 89)
@@ -1,6 +1,6 @@
-Date: 2015-06-09
+Date: 2015-06-13
Package: CHNOSZ
-Version: 1.0.5-1
+Version: 1.0.5-2
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/R/diagram.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -150,7 +150,6 @@
else if(as.residue & eout.is.aout) pv[[i]] <- pv[[i]] + eout$species$logact[i] / n.balance[i]
}
predominant <- which.pmax(pv)
- dim(predominant) <- dim(pv[[1]])
}
# a warning about that we can only show properties of the first species on a 2-D diagram
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/R/examples.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -30,7 +30,7 @@
}
demos <- function(which=c("sources", "NaCl", "density",
- "phosphate", "nucleobase", "orp", "diagram", "revisit", "findit",
+ "phosphate", "nucleobase", "ORP", "diagram", "revisit", "findit",
"CO2Ac", "nonideal", "ionize", "buffer", "yeastgfp", "mosaic",
"solubility", "wjd"), do.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
Modified: pkg/CHNOSZ/R/iprotein.R
===================================================================
--- pkg/CHNOSZ/R/iprotein.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/R/iprotein.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -171,7 +171,12 @@
# now we're ready to go
tp.new <- thermo$protein
if(!all(ipdup)) tp.new <- rbind(tp.new, aa[!ipdup, ])
- if(any(ipdup)) tp.new[ip[ipdup], ] <- aa[ipdup, ]
+ if(any(ipdup)) {
+ if(any(sapply(1:4, function(i){is.factor(aa[, i])})))
+ stop(paste("converting factors causes problems replacing protein data",
+ " data file should be read using e.g. aa <- read.csv(file, stringsAsFactors=FALSE)", sep="\n"))
+ tp.new[ip[ipdup], ] <- aa[ipdup, ]
+ }
rownames(tp.new) <- NULL
thermo$protein <- tp.new
assign("thermo", thermo, "CHNOSZ")
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/R/mosaic.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -58,12 +58,14 @@
if(blend) {
# calculate affinities using relative abundances of basis species
e <- equilibrate(A.bases)
+ # what is the total activity of the basis species?
+ loga.tot <- sum(10^unlist(e$loga.equil))
for(j in seq_along(affs)) {
for(i in seq_along(A.species$values)) {
# start with zero affinity
if(j==1) A.species$values[[i]][] <- 0
- # add affinity scaled by relative abundance of this basis species
- A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]
+ # add affinity scaled by __relative__ abundance of this basis species
+ A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]/loga.tot
}
}
} else {
Modified: pkg/CHNOSZ/R/util.list.R
===================================================================
--- pkg/CHNOSZ/R/util.list.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/R/util.list.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -3,6 +3,7 @@
which.pmax <- function (elts, na.rm = FALSE, pmin=FALSE) {
# adapted from R's pmax. elts is a list of numeric vectors
+ keepattr <- attributes(elts[[1]])
if(!is.numeric(elts[[1]])[1]) {
if(is.data.frame(elts[[1]])) elts[[1]] <- as.matrix(elts[[1]])
if(is.list(elts[[1]])) elts[[1]] <- elts[[1]][[1]]
@@ -11,29 +12,31 @@
mmm <- as.vector(elts[[1]])
which.mmm <- rep(1,length(elts[[1]]))
has.na <- FALSE
- for (i in 2:length(elts)) {
- if(!is.numeric(elts[[i]])[1]) {
- if(is.list(elts[[i]])) elts[[i]] <- elts[[i]][[1]]
- else elts[[i]] <- as.numeric(elts[[i]])
+ if(length(elts) > 1) {
+ for (i in 2:length(elts)) {
+ if(!is.numeric(elts[[i]])[1]) {
+ if(is.list(elts[[i]])) elts[[i]] <- elts[[i]][[1]]
+ else elts[[i]] <- as.numeric(elts[[i]])
+ }
+ work <- cbind(mmm, as.vector(elts[[i]]))
+ nas <- is.na(work)
+ if (has.na || (has.na <- any(nas))) {
+ work[, 1][nas[, 1]] <- work[, 2][nas[, 1]]
+ work[, 2][nas[, 2]] <- work[, 1][nas[, 2]]
+ }
+ if(pmin) change <- work[, 1] > work[, 2]
+ else change <- work[, 1] < work[, 2]
+ change <- change & !is.na(change)
+ work[, 1][change] <- work[, 2][change]
+ which.mmm[change] <- i
+ if (has.na && !na.rm) {
+ work[, 1][nas[, 1] | nas[, 2]] <- NA
+ which.mmm[nas[, 1] | nas[, 2]] <- NA
+ }
+ mmm <- work[, 1]
}
- work <- cbind(mmm, as.vector(elts[[i]]))
- nas <- is.na(work)
- if (has.na || (has.na <- any(nas))) {
- work[, 1][nas[, 1]] <- work[, 2][nas[, 1]]
- work[, 2][nas[, 2]] <- work[, 1][nas[, 2]]
- }
- if(pmin) change <- work[, 1] > work[, 2]
- else change <- work[, 1] < work[, 2]
- change <- change & !is.na(change)
- work[, 1][change] <- work[, 2][change]
- which.mmm[change] <- i
- if (has.na && !na.rm) {
- work[, 1][nas[, 1] | nas[, 2]] <- NA
- which.mmm[nas[, 1] | nas[, 2]] <- NA
- }
- mmm <- work[, 1]
}
- mostattributes(mmm) <- attributes(elts[[1]])
+ mostattributes(which.mmm) <- keepattr
which.mmm
}
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/demo/00Index 2015-06-13 17:27:56 UTC (rev 89)
@@ -3,7 +3,7 @@
density density of H2O, inverted from IAPWS-95 equations
phosphate phosphate speciation with pH, temperature and ionic strength
nucleobase relative stabilities of nucleobases and some amino acids
-orp oxidation-reduction potential of redox standards as a function of temperature
+ORP oxidation-reduction potential of redox standards as a function of temperature
diagram comparison of methods for calculating stability fields
revisit detailed example of usage of revisit()
findit detailed example of usage of findit()
Copied: pkg/CHNOSZ/demo/ORP.R (from rev 60, pkg/CHNOSZ/demo/orp.R)
===================================================================
--- pkg/CHNOSZ/demo/ORP.R (rev 0)
+++ pkg/CHNOSZ/demo/ORP.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -0,0 +1,160 @@
+# yell2010/orp.R 20100715 jmd
+# calculate the temperature dependence of
+# potentials vs. standard hydrogen electrode (SHE) of various electrodes (Ag/AgCl)
+# and ORP standards (ZoBell, Light's, (tri)iodide)
+# CHNOSZ provides functions subcrt() and convert()
+# used in this example
+#require(CHNOSZ)
+# Bard et al.'s fit to the potential
+# (Bard, Parson, Jordan, Standard Potentials In Aqueous Solution, 1985)
+AgAgCl.Bard <- function(T,high.T=TRUE) {
+ # we use the corrected high-T formula from wikipedia
+ if(high.T) return(0.23737 - 5.3783e-4 * T - 2.3728e-6 * T^2 - 2.2671e-9 * (T+273))
+ else return(0.23695 - 4.8564e-4 * T - 3.4205e-6 * T^2 - 5.869e-9 * (T+273))
+}
+# function to calculate the potential of Ag/AgCl vs. SHE
+# Ag(s) + Cl- = AgCl(s) + e-
+# logK = -pe - logaCl
+# pe = -logK - logmCl - loggamCl
+# ORP = RT/F * (logK - logmCl - loggamCl)
+AgAgCl <- function(T,mKCl=4) {
+ # mKCl is the molality of KCl in the electrolyte
+ # we take it as a first approximation to be equal to
+ # the molality of Cl- (and to the ionic strength)
+ logmCl <- log10(mKCl)
+ # get the logK for the reaction
+ logK <- subcrt(c("Ag","Cl-","AgCl","e-"),c(-1,-1,1,1),c("cr","aq","cr","aq"),T=T)$out$logK
+ # get the activity coefficient for Cl-
+ loggamCl <- subcrt("Cl-",T=T,IS=mKCl)$out[[1]]$loggam
+ # get the pe for the solution
+ pe <- -logK - logmCl - loggamCl
+ # convert that to Eh
+ Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ return(Eh)
+}
+ZoBell <- function(T) {
+ # doesn't work very well because we ignore the
+ # ferricyanide and ferrocyanide complexes
+ # Fe+2 = Fe+3 + e-
+ # logK = logaFe3 - logaFe2 - pe
+ # get the logK for the reaction
+ logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
+ # we use the recipe from standard methods (table 2580:II)
+ # 1.4080 g K4Fe(CN)6.3H2O -> 0.0033333 mol Fe+2
+ # 1.0975 g K3Fe(CN)6 -> 0.0033333 mol Fe+3
+ # 7.4555 g KCl -> 0.1 mol Cl-
+ logmFe2 <- logmFe3 <- log10(0.0033333)
+ # get the loggam for the iron species
+ loggamFe2 <- subcrt("Fe+2",T=T,IS=1)$out[[1]]$loggam
+ loggamFe3 <- subcrt("Fe+3",T=T,IS=1)$out[[1]]$loggam
+ # get the pe for the solution
+ pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
+ # convert to Eh
+ Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ return(Eh)
+}
+ZoBell.table <- function(T=NULL,which=NULL) {
+ # oxidation-reduction potential of ZoBell's solution
+ # from Standard Methods for Water and Wastewater or YSI
+ # (interpolated and/or extrapolated as necessary)
+ # standard methods (1997) table 2580:I
+ Eh.T.SMW <- 1:30
+ Eh.SMW <- c(0.481,0.479,0.476,0.474,0.472,0.47,0.468,0.465,0.463,0.461,
+ 0.459,0.457,0.454,0.452,0.45,0.448,0.446,0.443,0.441,0.439,0.437,
+ 0.435,0.432,0.43,0.428,0.426,0.424,0.421,0.419,0.417)
+ # from YSI (2005):
+ # Measuring ORP on YSI 6-Series Sondes: Tips, Cautions and Limitations
+ # NOTE: these values are vs. Ag/AgCl (4 M KCl)
+ Eh.T.YSI <- seq(-5,50,by=5)
+ Eh.YSI <- c(267.0,260.5,254.0,247.5,241.0,234.5,228.0,221.5,215.0,208.5,202.0,195.5)/1000
+ # spline function for each of the tables
+ SMW <- splinefun(Eh.T.SMW,Eh.SMW)
+ YSI <- splinefun(Eh.T.YSI,Eh.YSI)
+ # just one of the tables
+ Eh.fun <- get(which)
+ Eh.T <- get(paste("Eh.T",which,sep="."))
+ if(is.null(T)) T <- Eh.T
+ return(data.frame(T=T,Eh=Eh.fun(T)))
+}
+Light <- function(T) {
+ # this is going to look something like
+ # Fe+2 = Fe+3 + e-
+ # logK = logaFe3 - logaFe2 - pe
+ # get the logK for the reaction
+ logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
+ # we use the recipe from standard methods (table 2580:II)
+ # 39.21 g Fe(NH4)2(SO4)2(H2O)6 -> 0.1 mol Fe+2
+ # 48.22 g Fe(NH4)(SO4)2(H2O)12 -> 0.1 mol Fe+3
+ logmFe2 <- logmFe3 <- log10(0.1)
+ # get the loggam for the iron species
+ loggamFe2 <- subcrt("Fe+2",T=T,IS=0.2)$out[[1]]$loggam
+ loggamFe3 <- subcrt("Fe+3",T=T,IS=0.2)$out[[1]]$loggam
+ # get the pe for the solution
+ pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
+ # convert to Eh
+ Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ return(Eh)
+}
+Iodide.table <- function(T=NULL) {
+ # oxidation-reduction potential of Thermo's iodide solution
+ # from thermo instruction sheet 255218-001 (articlesFile_18739)
+ T.Iodide <- seq(0,50,5)
+ Eh.Iodide <- c(438,435,431,428,424,420,415,411,406,401,396)/1000
+ Iodide <- splinefun(T.Iodide,Eh.Iodide)
+ if(is.null(T)) T <- T.Iodide
+ return(data.frame(T=T,Eh=Iodide(T)))
+}
+Iodide <- function(T) {
+ # this is going to look something like
+ # 3I- = I3- + 2e-
+ # logK = -2pe + logaI3 - 3logaI
+ # get the logK for the reaction
+ logK <- subcrt(c("I-","I3-","e-"),c(-3,1,2),T=T)$out$logK
+ # could the activities be 0.1 M ... or something else?
+ logmI <- log10(2)
+ logmI3 <- log10(0.01)
+ # get the loggam for the iodine species
+ loggamI <- subcrt("I-",T=T,IS=0.2)$out[[1]]$loggam
+ loggamI3 <- subcrt("I3-",T=T,IS=0.2)$out[[1]]$loggam
+ # get the pe for the solution
+ pe <- ( -logK + logmI3 + loggamI3 - 3 * (logmI - loggamI) ) / 2
+ # convert to Eh
+ Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ return(Eh)
+}
+figure <- function() {
+ # make some figures
+ # the temperatures we're interested in
+ # in degrees C
+ T <- seq(0,100,5)
+ # temperature-Eh diagram for various electrodes
+ thermo.plot.new(ylim=c(0,0.8),xlim=c(0,100),
+ ylab=axis.label("Eh"),xlab=axis.label("T"))
+ # the Ag/AgCl electrode (Bard et al. fit)
+ points(T,AgAgCl.Bard(T),pch=0)
+ # the Ag/AgCl electrode (equilibrium calculations)
+ lines(T,AgAgCl(T))
+ # ZoBell's solution (SMW table 2580)
+ SMW <- ZoBell.table(which="SMW")
+ points(SMW$T,SMW$Eh,pch=1)
+ # ZoBell's solution (YSI tech report table)
+ YSI <- ZoBell.table(which="YSI")
+ # make these values referenced to SHE instead of Ag/AgCl
+ Eh.YSI <- YSI$Eh + AgAgCl(YSI$T)
+ points(YSI$T,Eh.YSI,pch=2)
+ # Light's solution (equilibrium values)
+ lines(T,Light(T))
+ # Light's solution (at 25 degrees only)
+ points(25,0.475 + 0.200,pch=3)
+ # Thermo's I-/I3- solution
+ Thermo <- Iodide.table()
+ points(Thermo$T,Thermo$Eh,pch=4)
+ # calculated I-/I3- values
+ lines(T,Iodide(T))
+ # add some labels
+ text(c(30,30,30,50),c(0.72,0.5,0.35,0.25),
+ c("Light","ZoBell","(Tri)Iodide","Ag/AgCl"))
+ title(main="Potentials vs standard hydrogen electrode (SHE)")
+}
+# finally, make the plot
+figure()
Modified: pkg/CHNOSZ/demo/density.R
===================================================================
--- pkg/CHNOSZ/demo/density.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/demo/density.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -7,8 +7,9 @@
blue <- "blue"
if(TPrange=="low") {
- T <- seq(300, 700, 5)
- P <- seq(1, 351, 5)
+ T <- seq(300, 800, 10)
+ #P <- seq(1, 600, 11.98)
+ P <- seq(0, 600, 12)
bias <- 1.68
} else {
# upper T,P limit for SUPCRT92: 2250 degC, 30000 bar
@@ -37,30 +38,31 @@
rho.mat <- matrix(rho.num, nrow=length(T), ncol=length(P))
# blueest for most dense, reddest for least dense
# bias is adjusted to white for the critical density
-ncol <- 100
+ncol <- 500
col <- colorRampPalette(c("red", "white", blue), bias=bias)(ncol)
# first make a background image (for debugging -
# will be visible only if some density calculations fail)
fill.mat <- matrix(0, nrow=length(T), ncol=length(P))
-image(T, P, fill.mat, col="black", xlab=axis.label("T", "K"), ylab=axis.label("P"))
+image(T, P, fill.mat, col="black", xlab=axis.label("T", "K"), ylab=axis.label("P"), useRaster=TRUE, yaxt="n")
+axis(2, at=c(1, seq(100, 600, 100)))
# now plot densities
-image(T, P, rho.mat, col=col, add=TRUE)
+image(T, P, rho.mat, col=col, add=TRUE, useRaster=TRUE)
# add a title and calculate saturation line
if(method=="IAPWS95") {
title(main=expression("Density of"~H[2]*O~"inverted from IAPWS-95 equations"))
- title(main=expression("Line calculated using auxiliary equations for saturation"), line=0.8)
+## title(main=expression("Line calculated using auxiliary equations for saturation"), line=0.8)
Psat <- convert(WP02.auxiliary("P.sigma", T), "bar")
} else if(method=="SUPCRT92") {
title(main=expression("Density of"~H[2]*O~"calculated using SUPCRT92"))
Psat <- water.SUPCRT92("Psat", T, "Psat")[,1]
}
-# plot saturation line
-lines(T, Psat, lwd=6)
-lines(T, Psat, lwd=3, col="gold")
+### plot saturation line
+##lines(T, Psat, lwd=6)
+##lines(T, Psat, lwd=3, col="gold")
# add a color key
if(TPrange=="low") {
- x <- c(333, 366, 370)
- y <- c(100, 300)
+ x <- c(355, 395, 402)
+ y <- c(170, 520)
} else if(TPrange=="high") {
x <- c(600, 780, 800)
y <- c(10000, 25000)
@@ -70,13 +72,13 @@
rect(x[1], ykey[1], x[2], rev(ykey)[1])
# label the extrema
rrange <- range(rho.num, na.rm=TRUE)
-text(x[3], ykey[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[1], 4)))), adj=0)
-text(x[3], rev(ykey)[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[2], 4)))), adj=0)
+text(x[3], ykey[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[1], 4)))), adj=0, col="white")
+text(x[3], rev(ykey)[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[2], 4)))), adj=0, col="white")
# label the critical density
rlevels <- seq(rrange[1], rrange[2], length.out=ncol+1)
rho.critical <- 322
icrit <- which.min(abs(rlevels-rho.critical))
-text(x[3], ykey[icrit], as.expression(substitute(x~kg/m^3~group("(", rho[c], ")"), list(x=rho.critical))), adj=0)
+text(x[3], ykey[icrit], as.expression(substitute(x~kg/m^3~group("(", rho[c], ")"), list(x=rho.critical))), adj=0, col="white")
#if(method=="IAPWS95") {
# # the saturation line is very accurate but not quite perfect;
Deleted: pkg/CHNOSZ/demo/orp.R
===================================================================
--- pkg/CHNOSZ/demo/orp.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/demo/orp.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -1,160 +0,0 @@
-# yell2010/orp.R 20100715 jmd
-# calculate the temperature dependence of
-# potentials vs. SHE of various electrodes (Ag/AgCl)
-# and ORP standards (ZoBell, Light's, (tri)iodide)
-# CHNOSZ provides functions subcrt() and convert()
-# used in this example
-#require(CHNOSZ)
-# Bard et al.'s fit to the potential
-# (Bard, Parson, Jordan, Standard Potentials In Aqueous Solution, 1985)
-AgAgCl.Bard <- function(T,high.T=TRUE) {
- # we use the corrected high-T formula from wikipedia
- if(high.T) return(0.23737 - 5.3783e-4 * T - 2.3728e-6 * T^2 - 2.2671e-9 * (T+273))
- else return(0.23695 - 4.8564e-4 * T - 3.4205e-6 * T^2 - 5.869e-9 * (T+273))
-}
-# function to calculate the potential of Ag/AgCl vs. SHE
-# Ag(s) + Cl- = AgCl(s) + e-
-# logK = -pe - logaCl
-# pe = -logK - logmCl - loggamCl
-# ORP = RT/F * (logK - logmCl - loggamCl)
-AgAgCl <- function(T,mKCl=4) {
- # mKCl is the molality of KCl in the electrolyte
- # we take it as a first approximation to be equal to
- # the molality of Cl- (and to the ionic strength)
- logmCl <- log10(mKCl)
- # get the logK for the reaction
- logK <- subcrt(c("Ag","Cl-","AgCl","e-"),c(-1,-1,1,1),c("cr","aq","cr","aq"),T=T)$out$logK
- # get the activity coefficient for Cl-
- loggamCl <- subcrt("Cl-",T=T,IS=mKCl)$out[[1]]$loggam
- # get the pe for the solution
- pe <- -logK - logmCl - loggamCl
- # convert that to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
- return(Eh)
-}
-ZoBell <- function(T) {
- # doesn't work very well because we ignore the
- # ferricyanide and ferrocyanide complexes
- # Fe+2 = Fe+3 + e-
- # logK = logaFe3 - logaFe2 - pe
- # get the logK for the reaction
- logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
- # we use the recipe from standard methods (table 2580:II)
- # 1.4080 g K4Fe(CN)6.3H2O -> 0.0033333 mol Fe+2
- # 1.0975 g K3Fe(CN)6 -> 0.0033333 mol Fe+3
- # 7.4555 g KCl -> 0.1 mol Cl-
- logmFe2 <- logmFe3 <- log10(0.0033333)
- # get the loggam for the iron species
- loggamFe2 <- subcrt("Fe+2",T=T,IS=1)$out[[1]]$loggam
- loggamFe3 <- subcrt("Fe+3",T=T,IS=1)$out[[1]]$loggam
- # get the pe for the solution
- pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
- return(Eh)
-}
-ZoBell.table <- function(T=NULL,which=NULL) {
- # oxidation-reduction potential of ZoBell's solution
- # from Standard Methods for Water and Wastewater or YSI
- # (interpolated and/or extrapolated as necessary)
- # standard methods (1997) table 2580:I
- Eh.T.SMW <- 1:30
- Eh.SMW <- c(0.481,0.479,0.476,0.474,0.472,0.47,0.468,0.465,0.463,0.461,
- 0.459,0.457,0.454,0.452,0.45,0.448,0.446,0.443,0.441,0.439,0.437,
- 0.435,0.432,0.43,0.428,0.426,0.424,0.421,0.419,0.417)
- # from YSI (2005):
- # Measuring ORP on YSI 6-Series Sondes: Tips, Cautions and Limitations
- # NOTE: these values are vs. Ag/AgCl (4 M KCl)
- Eh.T.YSI <- seq(-5,50,by=5)
- Eh.YSI <- c(267.0,260.5,254.0,247.5,241.0,234.5,228.0,221.5,215.0,208.5,202.0,195.5)/1000
- # spline function for each of the tables
- SMW <- splinefun(Eh.T.SMW,Eh.SMW)
- YSI <- splinefun(Eh.T.YSI,Eh.YSI)
- # just one of the tables
- Eh.fun <- get(which)
- Eh.T <- get(paste("Eh.T",which,sep="."))
- if(is.null(T)) T <- Eh.T
- return(data.frame(T=T,Eh=Eh.fun(T)))
-}
-Light <- function(T) {
- # this is going to look something like
- # Fe+2 = Fe+3 + e-
- # logK = logaFe3 - logaFe2 - pe
- # get the logK for the reaction
- logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
- # we use the recipe from standard methods (table 2580:II)
- # 39.21 g Fe(NH4)2(SO4)2(H2O)6 -> 0.1 mol Fe+2
- # 48.22 g Fe(NH4)(SO4)2(H2O)12 -> 0.1 mol Fe+3
- logmFe2 <- logmFe3 <- log10(0.1)
- # get the loggam for the iron species
- loggamFe2 <- subcrt("Fe+2",T=T,IS=0.2)$out[[1]]$loggam
- loggamFe3 <- subcrt("Fe+3",T=T,IS=0.2)$out[[1]]$loggam
- # get the pe for the solution
- pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
- return(Eh)
-}
-Iodide.table <- function(T=NULL) {
- # oxidation-reduction potential of Thermo's iodide solution
- # from thermo instruction sheet 255218-001 (articlesFile_18739)
- T.Iodide <- seq(0,50,5)
- Eh.Iodide <- c(438,435,431,428,424,420,415,411,406,401,396)/1000
- Iodide <- splinefun(T.Iodide,Eh.Iodide)
- if(is.null(T)) T <- T.Iodide
- return(data.frame(T=T,Eh=Iodide(T)))
-}
-Iodide <- function(T) {
- # this is going to look something like
- # 3I- = I3- + 2e-
- # logK = -2pe + logaI3 - 3logaI
- # get the logK for the reaction
- logK <- subcrt(c("I-","I3-","e-"),c(-3,1,2),T=T)$out$logK
- # could the activities be 0.1 M ... or something else?
- logmI <- log10(2)
- logmI3 <- log10(0.01)
- # get the loggam for the iodine species
- loggamI <- subcrt("I-",T=T,IS=0.2)$out[[1]]$loggam
- loggamI3 <- subcrt("I3-",T=T,IS=0.2)$out[[1]]$loggam
- # get the pe for the solution
- pe <- ( -logK + logmI3 + loggamI3 - 3 * (logmI - loggamI) ) / 2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
- return(Eh)
-}
-figure <- function() {
- # make some figures
- # the temperatures we're interested in
- # in degrees C
- T <- seq(0,100,5)
- # temperature-Eh diagram for various electrodes
- thermo.plot.new(ylim=c(0,0.8),xlim=c(0,100),
- ylab=axis.label("Eh"),xlab=axis.label("T"))
- # the Ag/AgCl electrode (Bard et al. fit)
- points(T,AgAgCl.Bard(T),pch=0)
- # the Ag/AgCl electrode (equilibrium calculations)
- lines(T,AgAgCl(T))
- # ZoBell's solution (SMW table 2580)
- SMW <- ZoBell.table(which="SMW")
- points(SMW$T,SMW$Eh,pch=1)
- # ZoBell's solution (YSI tech report table)
- YSI <- ZoBell.table(which="YSI")
- # make these values referenced to SHE instead of Ag/AgCl
- Eh.YSI <- YSI$Eh + AgAgCl(YSI$T)
- points(YSI$T,Eh.YSI,pch=2)
- # Light's solution (equilibrium values)
- lines(T,Light(T))
- # Light's solution (at 25 degrees only)
- points(25,0.475 + 0.200,pch=3)
- # Thermo's I-/I3- solution
- Thermo <- Iodide.table()
- points(Thermo$T,Thermo$Eh,pch=4)
- # calculated I-/I3- values
- lines(T,Iodide(T))
- # add some labels
- text(c(30,30,30,50),c(0.72,0.5,0.35,0.25),
- c("Light","ZoBell","(Tri)Iodide","Ag/AgCl"))
- title(main="Potentials vs SHE")
-}
-# finally, make the plot
-figure()
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/inst/NEWS 2015-06-13 17:27:56 UTC (rev 89)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.5-1 (2015-06-09)
+CHANGES IN CHNOSZ 1.0.5-2 (2015-06-13)
--------------------------------------
- Rewrite rho.IAPWS95() to be able to invert density from IAPWS-95
@@ -23,6 +23,13 @@
https://github.com/hadley/testthat/issues/129
https://github.com/hadley/testthat/issues/144
+- Fix bugs in which.pmax() that prevented proper assignment of
+ attributes in output, and functionality for lists of length 1.
+
+- mosaic() now correctly multiplies affinities by _relative_ abundances
+ of basis species when blend=TRUE. Thanks to Grayson Boyer for the bug
+ report that led to this fix.
+
CHANGES IN CHNOSZ 1.0.5 (2015-05-19)
------------------------------------
Modified: pkg/CHNOSZ/inst/tests/test-iprotein.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-iprotein.R 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/inst/tests/test-iprotein.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -5,7 +5,14 @@
test_that("basic searches and conversions work as expected", {
expect_equal(iprotein(c("LYSC_CHICK", "MYGPHYCA")), c(6, NA))
-
+ # factors causing problems again ...
+ f <- system.file("extdata/protein/DS11.csv", package="CHNOSZ")
+ aa <- read.csv(f)
+ # this adds the proteins
+ ip <- add.protein(aa)
+ # the replaces the proteins (with the same ones)
+ expect_error(ip <- add.protein(aa), "converting factors causes problems replacing protein data")
+ # ... should use read.csv(file, stringsAsFactors=FALSE)
})
test_that("errors and messages occur in some circumstances", {
Added: pkg/CHNOSZ/inst/tests/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-mosaic.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-mosaic.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -0,0 +1,17 @@
+context("mosaic")
+
+test_that("results are consistent with affinity()", {
+ basis(c("CO2", "H2O", "NH3", "O2"), c(0, 0, 0, 0))
+ species(c("alanine", "glycine"))
+ a <- affinity()
+ # this is a degenerate case because we only allow NH3 to swap for NH3, and CO2 for CO2;
+ # however it still exercises the affinity scaling and summing code
+ m1 <- mosaic("NH3", "CO2", blend=TRUE)
+ # this failed before we divided by loga.tot to get _relative_ abundances of basis species in mosaic.R
+ expect_equal(a$values, m1$A.species$values)
+ # the next call failed when which.pmax(), called by diagram(), choked on a list of length one
+ m2 <- mosaic("NH3", "CO2")
+ expect_equal(a$values, m2$A.species$values)
+})
+
+# TODO: test that basis specifications can be exchanged between bases and bases2 without altering output
Added: pkg/CHNOSZ/inst/tests/test-util.list.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-util.list.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-util.list.R 2015-06-13 17:27:56 UTC (rev 89)
@@ -0,0 +1,8 @@
+context("util.list")
+
+test_that("which.pmax() properly applies attributes, and also works for lists of length 1", {
+ testlist <- list(a=matrix(c(1,2,3,4)), b=matrix(c(4,3,2,1)))
+ testattr <- attributes(testlist[[1]])
+ expect_equal(attributes(which.pmax(testlist)), testattr)
+ expect_equal(as.numeric(which.pmax(testlist[1])), c(1, 1, 1, 1))
+})
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2015-06-09 10:21:09 UTC (rev 88)
+++ pkg/CHNOSZ/man/examples.Rd 2015-06-13 17:27:56 UTC (rev 89)
@@ -14,7 +14,7 @@
\usage{
examples(do.png = FALSE)
demos(which = c("sources", "NaCl", "density",
- "phosphate", "nucleobase", "orp", "diagram", "revisit", "findit",
+ "phosphate", "nucleobase", "ORP", "diagram", "revisit", "findit",
"CO2Ac", "nonideal", "ionize", "buffer", "yeastgfp", "mosaic",
"solubility", "wjd"), do.png=FALSE)
}
@@ -36,7 +36,7 @@
\code{density} \tab density of H2O, inverted from IAPWS-95 equations (\code{\link{rho.IAPWS95}}) \cr
\code{phosphate} \tab phosphate speciation with pH, temperature and ionic strength \cr
\code{nucleobase} \tab relative stabilities of nucleobases and some amino acids \cr
- \code{orp} \tab oxidation-reduction potential of redox standards as a function of temperature \cr
+ \code{ORP} \tab oxidation-reduction potential of redox standards as a function of temperature \cr
\code{diagram} \tab comparison of methods for calculating stability fields \cr
\code{revisit} \tab detailed example of usage of \code{\link{revisit}} \cr
\code{findit} \tab detailed example of usage of \code{\link{findit}} \cr
@@ -69,7 +69,7 @@
\examples{
\dontshow{data(thermo)}
-demos(c("orp", "NaCl"))
+demos(c("ORP", "NaCl"))
\dontshow{par(thermo$opar)}
\dontrun{
# use the following to run examples in all help topics
More information about the CHNOSZ-commits
mailing list