[CHNOSZ-commits] r203 - in pkg/CHNOSZ: . R demo inst man man/macros tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jun 4 19:28:46 CEST 2017
Author: jedick
Date: 2017-06-04 19:28:46 +0200 (Sun, 04 Jun 2017)
New Revision: 203
Added:
pkg/CHNOSZ/tests/testthat/test-water.lines.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/demo/mosaic.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/macros/macros.Rd
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/man/util.plot.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
Log:
water.lines() works for more combinations of variables
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/DESCRIPTION 2017-06-04 17:28:46 UTC (rev 203)
@@ -1,6 +1,6 @@
-Date: 2017-05-24
+Date: 2017-06-04
Package: CHNOSZ
-Version: 1.1.0-1
+Version: 1.1.0-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 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/R/diagram.R 2017-06-04 17:28:46 UTC (rev 203)
@@ -157,17 +157,26 @@
}
}
predominant <- which.pmax(pv)
- # for an Eh-pH or pe-pH diagram, clip plot to water stability region
- if(limit.water & eout$vars[1] == "pH" & eout$vars[2] %in% c("Eh", "pe")) {
- wl <- water.lines(xaxis=eout$vars[1], yaxis=eout$vars[2], T=eout$T, P=eout$P, xpoints=eout$vals[[1]], plot.it=FALSE)
- # for each x-point, find the y-values that are outside the water stability limits
- for(i in seq_along(wl$xpoints)) {
- ymin <- min(c(wl$y.oxidation[i], wl$y.reduction[i]))
- ymax <- max(c(wl$y.oxidation[i], wl$y.reduction[i]))
- # the actual calculation
- iNA <- eout$vals[[2]] < ymin | eout$vals[[2]] > ymax
- # assign NA to the predominance matrix
- predominant[i, iNA] <- NA
+ # clip plot to water stability region
+ if(limit.water) {
+ wl <- water.lines(eout, plot.it=FALSE)
+ # proceed if water.lines produced calculations for this plot
+ if(!identical(wl, NA)) {
+ # for each x-point, find the y-values that are outside the water stability limits
+ for(i in seq_along(wl$xpoints)) {
+ ymin <- min(c(wl$y.oxidation[i], wl$y.reduction[i]))
+ ymax <- max(c(wl$y.oxidation[i], wl$y.reduction[i]))
+ if(!wl$swapped) {
+ # the actual calculation
+ iNA <- eout$vals[[2]] < ymin | eout$vals[[2]] > ymax
+ # assign NA to the predominance matrix
+ predominant[i, iNA] <- NA
+ } else {
+ # as above, but x- and y-axes are swapped
+ iNA <- eout$vals[[1]] < ymin | eout$vals[[1]] > ymax
+ predominant[iNA, i] <- NA
+ }
+ }
}
}
}
Modified: pkg/CHNOSZ/R/util.plot.R
===================================================================
--- pkg/CHNOSZ/R/util.plot.R 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/R/util.plot.R 2017-06-04 17:28:46 UTC (rev 203)
@@ -65,67 +65,98 @@
par(opar)
}
-water.lines <- function(xaxis='pH', yaxis='Eh', T=298.15, P='Psat', which=c('oxidation','reduction'),
- logaH2O=0, lty=2, lwd=1, col=par('fg'), xpoints=NULL, O2state="gas", plot.it=TRUE) {
+water.lines <- function(eout, which=c('oxidation','reduction'),
+ lty=2, lwd=1, col=par('fg'), plot.it=TRUE) {
# draw water stability limits
# for Eh-pH, logfO2-pH, logfO2-T or Eh-T diagrams
# (i.e. redox variable is on the y axis)
+ # get axes, T, P, and xpoints from output of affinity() or equilibrate()
+ if(missing(eout)) stop("'eout' (the output of affinity(), equilibrate(), or diagram()) is missing")
+ # we only work on diagrams with 2 variables
+ if(length(eout$vars) != 2) return(NA)
+ # if needed, swap axes so T or P is on x-axis
+ swapped <- FALSE
+ if(eout$vars[2] %in% c("T", "P")) {
+ eout$vars <- rev(eout$vars)
+ eout$vals <- rev(eout$vals)
+ swapped <- TRUE
+ }
+ xaxis <- eout$vars[1]
+ yaxis <- eout$vars[2]
+ xpoints <- eout$vals[[1]]
+ # T and P are constants unless they are plotted on one of the axes
+ T <- eout$T
+ if(eout$vars[1]=="T") T <- envert(xpoints, "K")
+ P <- eout$P
+ if(eout$vars[1]=="P") P <- envert(xpoints, "bar")
+ # logaH2O is 0 unless given in eout$basis
+ iH2O <- match("H2O", rownames(eout$basis))
+ if(is.na(iH2O)) logaH2O <- 0 else logaH2O <- eout$basis$logact[iH2O]
+ # pH is 7 unless given in eout$basis or plotted on one of the axes
+ iHplus <- match("H+", rownames(eout$basis))
+ if(eout$vars[1]=="pH") pH <- xpoints
+ else if(!is.na(iHplus)) {
+ minuspH <- eout$basis$logact[iHplus]
+ # special treatment for non-numeric value (happens when a buffer is used, even for another basis species)
+ if(can.be.numeric(minuspH)) pH <- -as.numeric(minuspH) else pH <- NA
+ }
+ else pH <- 7
+ # O2state is gas unless given in eout$basis
+ iO2 <- match("O2", rownames(eout$basis))
+ if(is.na(iO2)) O2state <- "gas" else O2state <- eout$basis$state[iO2]
+ # H2state is gas unles given in eout$basis
+ iH2 <- match("H2", rownames(eout$basis))
+ if(is.na(iH2)) H2state <- "gas" else H2state <- eout$basis$state[iH2]
+ # where the calculated values will go
y.oxidation <- y.reduction <- NULL
- # if they are not provided, get the x points from the plot limits
- if(is.null(xpoints)) {
- pu <- par('usr')
- xlim <- pu[1:2]
- xpoints <- seq(xlim[1], xlim[2], length.out=100)
- }
- # note: Eh calculations are valid at a single T only
- if(xaxis=="O2" | (xaxis=='pH' & (yaxis=='Eh' | yaxis=='O2' | yaxis=="pe"))) {
+ if(xaxis %in% c("pH", "T", "P") & yaxis %in% c("Eh", "pe", "O2", "H2")) {
+ # Eh/pe/logfO2/logaO2/logfH2/logaH2 vs pH/T/P
if('reduction' %in% which) {
- logfH2 <- 0
- logK <- subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", O2state, "gas"), T=T, P=P, convert=FALSE)$out$logK
- # this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
- logfO2 <- 2 * logK - logfH2 + 2 * logaH2O
- #if(xaxis=='O2') abline(v=logfO2,lty=lty,lwd=lwd,col=col)
- #if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
- if(yaxis=="Eh") y.reduction <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
- else if(yaxis=="pe") y.reduction <- convert(convert(logfO2, 'E0', T=T, P=P, pH=xpoints), "pe", T=T)
+ logfH2 <- logaH2O # usually 0
+ if(yaxis=="H2") {
+ logK <- subcrt(c("H2", "H2"), c(-1, 1), c("gas", H2state), T=T, P=P, convert=FALSE)$out$logK
+ # this is logfH2 if H2state=="gas", or logaH2 if H2state=="aq"
+ logfH2 <- logfH2 + logK
+ y.reduction <- rep(logfH2, length.out=length(xpoints))
+ } else {
+ logK <- subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", O2state, "gas"), T=T, P=P, convert=FALSE)$out$logK
+ # this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
+ logfO2 <- 2 * (logK - logfH2 + logaH2O)
+ if(yaxis=="O2") y.reduction <- rep(logfO2, length.out=length(xpoints))
+ else if(yaxis=="Eh") y.reduction <- convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O)
+ else if(yaxis=="pe") y.reduction <- convert(convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O), "pe", T=T)
+ }
}
if('oxidation' %in% which) {
- logfO2 <- 0
- logK <- subcrt(c("O2", "O2"), c(-1, 1), c("gas", O2state), T=T, P=P, convert=FALSE)$out$logK
- # this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
- logfO2 <- logfO2 + logK
- #if(xaxis=='O2') abline(v=logfO2,lty=lty,lwd=lwd,col=col)
- #if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
- if(yaxis=="Eh") y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
- else if(yaxis=="pe") y.oxidation <- convert(convert(logfO2, 'E0', T=T, P=P, pH=xpoints), "pe", T=T)
+ logfO2 <- logaH2O # usually 0
+ if(yaxis=="H2") {
+ logK <- subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", "gas", H2state), T=T, P=P, convert=FALSE)$out$logK
+ # this is logfH2 if H2state=="gas", or logaH2 if H2state=="aq"
+ logfH2 <- logK - 0.5*logfO2 + logaH2O
+ y.oxidation <- rep(logfH2, length.out=length(xpoints))
+ } else {
+ logK <- subcrt(c("O2", "O2"), c(-1, 1), c("gas", O2state), T=T, P=P, convert=FALSE)$out$logK
+ # this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
+ logfO2 <- logfO2 + logK
+ if(yaxis=="O2") y.oxidation <- rep(logfO2, length.out=length(xpoints))
+ else if(yaxis=="Eh") y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O)
+ else if(yaxis=="pe") y.oxidation <- convert(convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O), "pe", T=T)
+ }
}
- } else if(xaxis %in% c('T','P') & yaxis %in% c('Eh','O2') ) {
- #if(xaxis=='T') if(is.null(xpoints)) xpoints <- T
- # 20090212 get T values from plot limits
- # TODO: make this work for T on y-axis too
- if(xaxis=='T' & missing(T)) T <- envert(xpoints, "K")
- if(xaxis=='P') if(missing(xpoints)) xpoints <- P
- if('oxidation' %in% which) {
- logfO2 <- rep(0,length(xpoints))
- if(yaxis=='Eh') y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
- else y.oxidation <- logfO2
- }
- if('reduction' %in% which) {
- logfH2 <- 0
- logK <- subcrt(c('H2O','oxygen','hydrogen'),c(-1,0.5,1),T=T,P=P,convert=FALSE)$out$logK
- logfO2 <- 2 * logK - logfH2 + 2 * logaH2O
- if(yaxis=='Eh') y.reduction <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
- else y.reduction <- logfO2
- }
- }
- if(yaxis=="Eh") y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
+ } else return(NA)
# now plot the lines
if(plot.it) {
- lines(xpoints, y.oxidation, lty=lty, lwd=lwd, col=col)
- lines(xpoints, y.reduction, lty=lty, lwd=lwd, col=col)
+ if(swapped) {
+ # xpoints above is really the ypoints
+ lines(y.oxidation, xpoints, lty=lty, lwd=lwd, col=col)
+ lines(y.reduction, xpoints, lty=lty, lwd=lwd, col=col)
+ } else {
+ lines(xpoints, y.oxidation, lty=lty, lwd=lwd, col=col)
+ lines(xpoints, y.reduction, lty=lty, lwd=lwd, col=col)
+ }
}
# return the values
- return(invisible(list(xpoints=xpoints, y.oxidation=y.oxidation, y.reduction=y.reduction)))
+ return(invisible(list(xpoints=xpoints, y.oxidation=y.oxidation, y.reduction=y.reduction, swapped=swapped)))
}
mtitle <- function(main, line=0, ...) {
Modified: pkg/CHNOSZ/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/demo/copper.R 2017-06-04 17:28:46 UTC (rev 203)
@@ -58,7 +58,7 @@
text(d$lx, -0.5, Gly, col="darkblue")
# add water lines and title
-water.lines()
+water.lines(d)
mtitle(expression("Copper-water-glycine at 25"~degree*"C and 1 bar",
"After Aksu and Doyle, 2001 (Fig. 2b)"), line=0.5)
Modified: pkg/CHNOSZ/demo/mosaic.R
===================================================================
--- pkg/CHNOSZ/demo/mosaic.R 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/demo/mosaic.R 2017-06-04 17:28:46 UTC (rev 203)
@@ -25,7 +25,7 @@
m1 <- mosaic(bases, bases2, TRUE, pH=pH, Eh=Eh, T=T)
# make a diagram and add water stability lines
diagram(m1$A.species, lwd=2)
-water.lines("pH", "Eh", T=convert(T, "K"), col="seagreen", lwd=1.5)
+water.lines(m1$A.species, col="seagreen", lwd=1.5)
# show lines for Fe(aq) = 10^-4 M
species(c("Fe+2", "Fe+3"), -4)
m2 <- mosaic(bases, bases2, TRUE, pH=pH, Eh=Eh, T=T)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/inst/NEWS 2017-06-04 17:28:46 UTC (rev 203)
@@ -1,6 +1,10 @@
-CHANGES IN CHNOSZ 1.1.0-1 (2017-05-24)
+CHANGES IN CHNOSZ 1.1.0-2 (2017-06-04)
--------------------------------------
+- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
+ logfH2, or logaH2 vs pH, T, or P. It is possible to have T or P on
+ either the x- or y-axis.
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/man/diagram.Rd 2017-06-04 17:28:46 UTC (rev 203)
@@ -218,7 +218,7 @@
d <- diagram(a, main="Fe-S-O-H, after Majzlan et al., 2006")
# the first four species show up in order near pe=15
stopifnot(all.equal(unique(d$predominant[, 183]), 1:4))
-water.lines(yaxis="pe", lwd=2)
+water.lines(d, lwd=2)
text(3, 22, describe.basis(thermo$basis[2:3,], digits=2, oneline=TRUE))
text(3, 21, describe.property(c("T", "P"), c(25, 1), oneline=TRUE))
Modified: pkg/CHNOSZ/man/macros/macros.Rd
===================================================================
--- pkg/CHNOSZ/man/macros/macros.Rd 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/man/macros/macros.Rd 2017-06-04 17:28:46 UTC (rev 203)
@@ -1,3 +1,4 @@
+% chemical formulas
\newcommand{\CO2}{\ifelse{latex}{\eqn{\mathrm{CO_{2}}}}{\ifelse{html}{\out{CO<sub>2</sub>}}{CO2}}}
\newcommand{\NH3}{\ifelse{latex}{\eqn{\mathrm{NH_{3}}}}{\ifelse{html}{\out{NH<sub>3</sub>}}{NH3}}}
\newcommand{\H2S}{\ifelse{latex}{\eqn{\mathrm{H_{2}S}}}{\ifelse{html}{\out{H<sub>2</sub>S}}{H2S}}}
@@ -3,9 +4,7 @@
\newcommand{\H2O}{\ifelse{latex}{\eqn{\mathrm{H_{2}O}}}{\ifelse{html}{\out{H<sub>2</sub>O}}{H2O}}}
\newcommand{\O2}{\ifelse{latex}{\eqn{\mathrm{O_{2}}}}{\ifelse{html}{\out{O<sub>2</sub>}}{O2}}}
-
\newcommand{\Hplus}{\ifelse{latex}{\eqn{\mathrm{H^{+}}}}{\ifelse{html}{\out{H<sup>+</sup>}}{H+}}}
\newcommand{\eminus}{\ifelse{latex}{\eqn{e^{-}}}{\ifelse{html}{\out{<i>e</i><sup>-</sup>}}{e-}}}
\newcommand{\Mgplus2}{\ifelse{latex}{\eqn{\mathrm{Mg^{+2}}}}{\ifelse{html}{\out{Mg<sup>+2</sup>}}{Mg+2}}}
-
\newcommand{\H3PO4}{\ifelse{latex}{\eqn{\mathrm{H_{3}PO_{4}}}}{\ifelse{html}{\out{H<sub>3</sub>PO<sub>4</sub>}}{H3PO4}}}
\newcommand{\Fe2O3}{\ifelse{latex}{\eqn{\mathrm{Fe_{2}O_{3}}}}{\ifelse{html}{\out{Fe<sub>2</sub>O<sub>3</sub>}}{Fe2O3}}}
@@ -14,13 +13,13 @@
\newcommand{\NH4plus}{\ifelse{latex}{\eqn{\mathrm{NH_{4}^{+}}}}{\ifelse{html}{\out{NH<sub>4</sub><sup>+</sup>}}{NH4+}}}
\newcommand{\H2}{\ifelse{latex}{\eqn{\mathrm{H_{2}}}}{\ifelse{html}{\out{H<sub>2</sub>}}{H2}}}
-
-
-
+% other common variables
+\newcommand{\T}{\ifelse{latex}{\eqn{T}}{\ifelse{html}{\out{<I>T</I>}}{T}}}
+\newcommand{\P}{\ifelse{latex}{\eqn{P}}{\ifelse{html}{\out{<I>P</I>}}{P}}}
+\newcommand{\logfO2}{\ifelse{latex}{\eqn{\log f_{\mathrm{O_{2}}}}}{\ifelse{html}{\out{log<i>f</i><sub>O<sub>2</sub></sub>}}{logfO2}}}
+\newcommand{\logfH2}{\ifelse{latex}{\eqn{\log f_{\mathrm{H_{2}}}}}{\ifelse{html}{\out{log<i>f</i><sub>H<sub>2</sub></sub>}}{logfH2}}}
+\newcommand{\logaH2O}{\ifelse{latex}{\eqn{\log a_{\mathrm{H_{2}O}}}}{\ifelse{html}{\out{log<i>a</i><sub>H<sub>2</sub>O</sub>}}{logaH2O}}}
\newcommand{\ZC}{\ifelse{latex}{\eqn{Z_\mathrm{C}}}{\ifelse{html}{\out{<I>Z</I><sub>C</sub>}}{ZC}}}
% use \nH2O{̄} to call this macro (the html code can't be defined in the macro,
% which interprets '#' followed by a number as a placeholder for an argument)
\newcommand{\nH2O}{\ifelse{latex}{\eqn{\bar{n}_{\mathrm{H_{2}O}}}}{\ifelse{html}{\out{<i>n</i>#1<sub>H<sub>2</sub>O</sub>}}{nH2O}}}
-
-\newcommand{\logfO2}{\ifelse{latex}{\eqn{\log f_{\mathrm{O_{2}}}}}{\ifelse{html}{\out{log<i>f</i><sub>O<sub>2</sub></sub>}}{logfO2}}}
-\newcommand{\logaH2O}{\ifelse{latex}{\eqn{\log a_{\mathrm{H_{2}O}}}}{\ifelse{html}{\out{log<i>a</i><sub>H<sub>2</sub>O</sub>}}{logaH2O}}}
Modified: pkg/CHNOSZ/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/man/mosaic.Rd 2017-06-04 17:28:46 UTC (rev 203)
@@ -63,8 +63,8 @@
# calculate affinities using the predominant basis species
m1 <- mosaic(bases, pH=pH, Eh=Eh, T=T)
# make a diagram and add water stability lines
-diagram(m1$A.species, lwd=2)
-water.lines("pH", "Eh", T=convert(T, "K"), col="seagreen", lwd=1.5)
+d <- diagram(m1$A.species, lwd=2)
+water.lines(d, col="seagreen", lwd=1.5)
# show lines for Fe(aq) = 10^-4 M
species(c("Fe+2", "Fe+3"), -4)
m2 <- mosaic(bases, pH=pH, Eh=Eh, T=T)
Modified: pkg/CHNOSZ/man/util.plot.Rd
===================================================================
--- pkg/CHNOSZ/man/util.plot.Rd 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/man/util.plot.Rd 2017-06-04 17:28:46 UTC (rev 203)
@@ -14,7 +14,6 @@
Initialize a new plot window using preset parameters, add an axis or title to a plot, generate labels for axes and subplots, add stability lines for water, get colors for a set of numeric values.
}
-
\usage{
thermo.plot.new(xlim, ylim, xlab, ylab, cex = par("cex"),
mar = NULL, lwd = par("lwd"), side = c(1,2,3,4),
@@ -26,9 +25,8 @@
usrfig()
label.figure(x, xfrac = 0.05, yfrac = 0.95, paren = FALSE,
italic = FALSE, ...)
- water.lines(xaxis = "pH", yaxis = "Eh", T = 298.15, P = "Psat",
- which = c("oxidation","reduction"), logaH2O = 0, lty = 2, lwd=1,
- col = par("fg"), xpoints = NULL, O2state="gas", plot.it = TRUE)
+ water.lines(eout, which = c("oxidation","reduction"),
+ lty = 2, lwd=1, col = par("fg"), plot.it = TRUE)
mtitle(main, line=0, ...)
ZC.col(z)
}
@@ -51,20 +49,14 @@
\item{las}{numeric, style for axis labels}
\item{xline}{numeric, margin line on which to plot \eqn{x}{x}-axis name}
\item{...}{further arguments passed to \code{par} or \code{mtext}}
- \item{T}{numeric, temperature (K)}
\item{x}{character, label to place on plot}
\item{xfrac}{numeric, fractional location on \eqn{x}{x}-axis for placement of label}
\item{yfrac}{numeric, fractional location on \eqn{y}{y}-axis for placement of label}
\item{paren}{logical, add parentheses around label text?}
\item{italic}{logical, italicize label text?}
- \item{xaxis}{character, description of \eqn{x}{x}-axis}
- \item{yaxis}{character, description of \eqn{y}{y}-axis}
- \item{P}{numeric, pressure (bar)}
+ \item{eout}{data frame, output of \code{\link{affinity}}, \code{\link{equilibrate}}, or \code{\link{diagram}}}
\item{which}{character, which of oxidation/reduction lines to plot}
- \item{logaH2O}{numeric, logarithm of the activity of \eqn{\mathrm{H_2O}}{H2O}}
\item{lty}{numeric, line type}
- \item{xpoints}{numeric, points to plot on \eqn{x}{x} axis}
- \item{O2state}{character, state of O2}
\item{plot.it}{logical, plot the lines?}
\item{main}{character, text for plot title}
\item{line}{numeric, margin line to place title}
@@ -75,12 +67,12 @@
\code{thermo.plot.new} sets parameters for a new plot, creates a new plot using \code{\link{plot.new}}, and adds axes and major and minor tick marks to the plot. Plot parameters (see \code{\link{par}}) including \code{cex}, \code{mar}, \code{lwd}, \code{mgp} and \code{axs} can be given, as well as a numeric vector in \code{side} identifying which sides of the plot receive tick marks. \code{yline}, if present, denotes the margin line (default \code{\link{par}('mgp')[1]}) where the y-axis name is plotted.
-\code{water.lines} plots lines representing the oxidation and reduction stability limits of water on \code{yaxis}-\code{xaxis} diagrams, where \code{yaxis} can be \samp{Eh} or \samp{O2}, and \code{xaxis} can be \samp{pH} or \samp{T}.
+\code{water.lines} plots lines representing the oxidation and reduction stability limits of water on Eh/pe/\logfO2/\logfH2 vs pH/\T/\P diagrams.
+The x- and y-variables and their ranges are taken from \code{eout}.
+Values of \T, \P, pH, and \logaH2O, not corresponding to either axis, are also taken from \code{eout}.
\code{which} controls which lines are drawn (\samp{oxidation}, \samp{reduction}, or both (the default)).
-\code{logaH2O} denotes the logarithm of the activity of water.
-With \code{O2state} set to \samp{gas} (the default), the logarithm of oxygen fugacity is plotted.
-Change this to \samp{aq} to plot the logarithm of oxygen activity (do not change it if plotting \samp{Eh}).
-\code{xpoints} is an optional list of points on the x axis to which to restrict the plotting (default of \code{NULL} refers to the axis limits).
+The value of \code{swapped} in the output reflects whether pH, \T, or \P is on the x-axis (TRUE) or y-axis (FALSE).
+\code{NA} is returned for any diagram for variables that can not be processed (including diagrams with more than 2 variables).
\code{label.plot} and \code{label.figure} add identifying text within the plot region and figure region.
The value given for \code{x} is made into a label, optionally italicized and with parentheses (like \ifelse{latex}{\eqn{(a)}}{\ifelse{html}{\out{(<i>a</i>)}}{(a)}}).
Added: pkg/CHNOSZ/tests/testthat/test-water.lines.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-water.lines.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-water.lines.R 2017-06-04 17:28:46 UTC (rev 203)
@@ -0,0 +1,73 @@
+# testing revised water.lines 2017-06-04
+context("water.lines")
+
+# change to TRUE to show plots (for additional testing)
+plot.it <- FALSE
+# change to FALSE to make base plots (for additional testing)
+limit.water <- TRUE
+
+# function to count the number of species on a diagram (not including NA fields)
+nspecies <- function(a) length(na.omit(unique(as.numeric(diagram(a, limit.water=limit.water, plot.it=plot.it)$predominant))))
+
+res <- 25
+
+test_that("water.lines() masks H2 and O2 fields on diagrams with pH as x-axis", {
+ if(plot.it) par(mfrow=c(3, 6))
+ Ts <- c(25, 100, 25)
+ logaH2O <- c(0, 0, -5)
+ for(i in 1:3) {
+ T <- Ts[i]
+ basis(c("H2O", "H+", "e-"), c(logaH2O[i], -7, -7))
+ species(c("H2O", "O2", "H2"), c("liq", "gas", "gas"))
+ # Eh-pH
+ n1 <- nspecies(affinity(pH=c(-2, 16, res), Eh=c(-1, 1.5, res), T=T))
+ # pe-pH
+ n2 <- nspecies(affinity(pH=c(-2, 16, res), pe=c(-20, 25, res), T=T))
+ # logfO2-pH
+ swap.basis("e-", "oxygen")
+ n3 <- nspecies(affinity(pH=c(-2, 16, res), O2=c(-90, 10, res), T=T))
+ # logaO2-pH
+ swap.basis("oxygen", "O2")
+ n4 <- nspecies(affinity(pH=c(-2, 16, res), O2=c(-90, 10, res), T=T))
+ # logfH2-pH
+ swap.basis("O2", "hydrogen")
+ n5 <- nspecies(affinity(pH=c(-2, 16, res), H2=c(-50, 10, res), T=T))
+ # logaH2-pH
+ swap.basis("hydrogen", "H2")
+ n6 <- nspecies(affinity(pH=c(-2, 16, res), H2=c(-50, 10, res), T=T))
+ expect_equal(c(n1, n2, n3, n4, n5, n6), c(1, 1, 1, 1, 1, 1))
+ }
+})
+
+if(plot.it) par(mfrow=c(3, 6))
+basis(c("H2O", "H+", "e-"), c(0, -7, -7))
+species(c("H2O", "O2", "H2"), c("liq", "gas", "gas"))
+test_that("water.lines() masks H2 and O2 fields on diagrams with T as x-axis", {
+ n1 <- nspecies(affinity(T=c(0, 200, res), Eh=c(-1, 1.5, res))) # Eh-T
+ n2 <- nspecies(affinity(T=c(0, 200, res), pe=c(-20, 25, res))) # pe-T
+ swap.basis("e-", "oxygen"); n3 <- nspecies(affinity(T=c(0, 200, res), O2=c(-90, 10, res))) # logfO2-T
+ swap.basis("oxygen", "O2"); n4 <- nspecies(affinity(T=c(0, 200, res), O2=c(-90, 10, res))) # logaO2-T
+ swap.basis("O2", "hydrogen"); n5 <- nspecies(affinity(T=c(0, 200, res), H2=c(-50, 10, res))) # logfH2-T
+ swap.basis("hydrogen", "H2"); n6 <- nspecies(affinity(T=c(0, 200, res), H2=c(-50, 10, res))) # logaH2-T
+ expect_equal(c(n1, n2, n3, n4, n5, n6), c(1, 1, 1, 1, 1, 1))
+})
+test_that("water.lines() masks H2 and O2 fields on diagrams with P as x-axis", {
+ swap.basis("H2", "e-")
+ n1 <- nspecies(affinity(P=c(1, 5000, res), Eh=c(-1, 1.5, res))) # Eh-P
+ n2 <- nspecies(affinity(P=c(1, 5000, res), pe=c(-20, 25, res))) # pe-P
+ swap.basis("e-", "oxygen"); n3 <- nspecies(affinity(P=c(1, 5000, res), O2=c(-90, 10, res))) # logfO2-P
+ swap.basis("oxygen", "O2"); n4 <- nspecies(affinity(P=c(1, 5000, res), O2=c(-90, 10, res))) # logaO2-P
+ swap.basis("O2", "hydrogen"); n5 <- nspecies(affinity(P=c(1, 5000, res), H2=c(-50, 10, res))) # logfH2-P
+ swap.basis("hydrogen", "H2"); n6 <- nspecies(affinity(P=c(1, 5000, res), H2=c(-50, 10, res))) # logaH2-P
+ expect_equal(c(n1, n2, n3, n4, n5, n6), c(1, 1, 1, 1, 1, 1))
+})
+test_that("water.lines() masks H2 and O2 fields on diagrams with T as y-axis", {
+ swap.basis("H2", "e-")
+ n1 <- nspecies(affinity(Eh=c(-1, 1.5, res), T=c(0, 200, res))) # T-Eh
+ n2 <- nspecies(affinity(pe=c(-20, 25, res), T=c(0, 200, res))) # T-pe
+ swap.basis("e-", "oxygen"); n3 <- nspecies(affinity(O2=c(-90, 10, res), T=c(0, 200, res))) # T-logfO2
+ swap.basis("oxygen", "O2"); n4 <- nspecies(affinity(O2=c(-90, 10, res), T=c(0, 200, res))) # T-logaO2
+ swap.basis("O2", "hydrogen"); n5 <- nspecies(affinity(H2=c(-50, 10, res), T=c(0, 200, res))) # T-logfH2
+ swap.basis("hydrogen", "H2"); n6 <- nspecies(affinity(H2=c(-50, 10, res), T=c(0, 200, res))) # T-logaH2
+ expect_equal(c(n1, n2, n3, n4, n5, n6), c(1, 1, 1, 1, 1, 1))
+})
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-06-04 17:28:46 UTC (rev 203)
@@ -547,7 +547,7 @@
```{r EhpH_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="Aqueous sulfur species at 25 °C.", pngquant=pngquant, timeit=timeit}
a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
diagram(a, fill = "heat")
-water.lines()
+water.lines(a)
```
Now we can calculate the affinities on an Eh-pH grid:
@@ -631,7 +631,7 @@
diagram(m1$A.species, lwd = 2, fill = NA, limit.water = FALSE)
diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 2,
limit.water = FALSE)
-water.lines("pH", "Eh", T = convert(T, "K"), col = "red", lwd = 2, lty = 2)
+water.lines(m1$A.species, col = "red", lwd = 2, lty = 2)
```
The argument `blend = TRUE` is used to combine the diagrams according to the equilibrium activities of the basis species by themselves ([see below](#equilibration)).
Modified: pkg/CHNOSZ/vignettes/equilibrium.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.Rnw 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/vignettes/equilibrium.Rnw 2017-06-04 17:28:46 UTC (rev 203)
@@ -1,4 +1,4 @@
-%% LyX 2.2.2 created this file. For more info, see http://www.lyx.org/.
+%% LyX 2.2.3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,round]{article}
\usepackage{mathpazo}
@@ -609,7 +609,7 @@
iCu <- which(d$predominant == 1, arr.ind=TRUE)
text(a$vals[[1]][max(iCu[, 1])] - 0.03, a$vals[[2]][min(iCu[, 2])] + 0.2, adj=1, loga)
}
-water.lines()
+water.lines(d)
@
@@ -1018,11 +1018,11 @@
layout(matrix(1:2), heights=c(1, 2))
e <- equilibrate(a)
diagram(e, ylim=c(-2.8, -1.6), names=organisms)
-water.lines(xaxis="O2")
+water.lines(e)
title(main="Equilibrium activities of proteins, normalize = FALSE")
e <- equilibrate(a, normalize=TRUE)
diagram(e, ylim=c(-4, -2), names=organisms)
-water.lines(xaxis="O2")
+water.lines(e)
title(main="Equilibrium activities of proteins, normalize = TRUE")
@
@@ -1048,7 +1048,7 @@
for(normalize in c(FALSE, TRUE)) {
e <- equilibrate(a, loga.balance=-2, normalize=normalize)
diagram(e, ylim=c(-30, 0), legend.x=NULL, col=col, lty=lty)
- water.lines(xaxis="O2", T=convert(325, "K"))
+ water.lines(e)
title(main=paste("Aqueous sulfur speciation, normalize =", normalize))
}
par(mar=c(0, 0, 0, 0))
Modified: pkg/CHNOSZ/vignettes/equilibrium.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.lyx 2017-05-24 16:06:11 UTC (rev 202)
+++ pkg/CHNOSZ/vignettes/equilibrium.lyx 2017-06-04 17:28:46 UTC (rev 203)
@@ -2685,7 +2685,7 @@
\begin_layout Plain Layout
-water.lines()
+water.lines(d)
\end_layout
\begin_layout Plain Layout
@@ -4233,7 +4233,7 @@
\begin_layout Plain Layout
-water.lines(xaxis="O2")
+water.lines(e)
\end_layout
\begin_layout Plain Layout
@@ -4253,7 +4253,7 @@
\begin_layout Plain Layout
-water.lines(xaxis="O2")
+water.lines(e)
\end_layout
\begin_layout Plain Layout
@@ -4361,7 +4361,7 @@
\begin_layout Plain Layout
- water.lines(xaxis="O2", T=convert(325, "K"))
+ water.lines(e)
\end_layout
\begin_layout Plain Layout
More information about the CHNOSZ-commits
mailing list