[CHNOSZ-commits] r190 - in pkg/CHNOSZ: . R demo inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 28 10:30:52 CEST 2017
Author: jedick
Date: 2017-04-28 10:30:51 +0200 (Fri, 28 Apr 2017)
New Revision: 190
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/demo/mosaic.R
pkg/CHNOSZ/demo/sources.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/util.plot.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
diagram(): gray areas beyond water stability limits
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/DESCRIPTION 2017-04-28 08:30:51 UTC (rev 190)
@@ -1,6 +1,6 @@
-Date: 2017-04-27
+Date: 2017-04-28
Package: CHNOSZ
-Version: 1.0.8-78
+Version: 1.0.8-79
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-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/R/diagram.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -20,7 +20,8 @@
# line styles
lty=NULL, lwd=par("lwd"), dotted=NULL,
# colors
- col=par("col"), col.names=par("col"), fill=NULL, col.NA="black",
+ col=par("col"), col.names=par("col"), fill=NULL,
+ fill.NA="slategray1", limit.water=TRUE,
# labels
names=NULL, main=NULL, legend.x=NA, format.names=TRUE, adj=0.5, dy=0,
# plotting controls
@@ -33,7 +34,7 @@
if(!"sout" %in% names(eout)) stop("'eout' does not look like output from equil() or affinity()")
## 'what' can be:
- # loga.equil - equilibrium activities of species of interest (eout)
+ # loga.equil - equilibrium activities of species of interest (eout)
# basis species - equilibrium activity of a basis species (aout)
# missing - property from affinity() or predominances of species (aout)
eout.is.aout <- FALSE
@@ -156,6 +157,19 @@
}
}
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
+ }
+ }
}
# a warning about that we can only show properties of the first species on a 2-D diagram
@@ -316,9 +330,9 @@
for(i in 1:nrow(zs)) zs[i,] <- out[nrow(zs)+1-i,]
zs <- t(zs)
breaks <- c(-1, 0, 1:nspecies) + 0.5
- # use col.NA for NA values
+ # use fill.NA for NA values
zs[is.na(zs)] <- 0
- image(x=xs, y=ys, z=zs, col=c(col.NA, fill), add=TRUE, breaks=breaks, useRaster=TRUE)
+ image(x=xs, y=ys, z=zs, col=c(fill.NA, fill), add=TRUE, breaks=breaks, useRaster=TRUE)
}
## curve plot function
# 20091116 replaced plot.curve with plot.line; different
@@ -405,7 +419,7 @@
ys <- rev(ys)
}
# the categories (species/groups/etc) on the plot
- zvals <- unique(as.vector(predominant))
+ zvals <- na.omit(unique(as.vector(predominant)))
# take each possible pair
for(i in 1:(length(zvals)-1)) {
for(j in (i+1):length(zvals)) {
@@ -460,6 +474,8 @@
else if(isTRUE(fill[1]=="rainbow")) fill <- rainbow(ngroups)
else if(isTRUE(fill[1] %in% c("heat", "terrain", "topo", "cm"))) fill <- get(paste0(fill[1], ".colors"))(ngroups)
fill <- rep(fill, length.out=ngroups)
+ # modify the default for fill.NA
+ if(add & missing(fill.NA)) fill.NA <- "transparent"
# the x and y values
xs <- eout$vals[[1]]
ys <- eout$vals[[2]]
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/R/mosaic.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -77,6 +77,8 @@
# merge affinities using the second, third, ... basis species
for(j in tail(seq_along(affs), -1)) {
is.predominant <- d$predominant==j
+ # diagram() produces NA beyond water limits on Eh-pH diagrams (but we can't use NA for indexing, below)
+ is.predominant[is.na(is.predominant)] <- FALSE
for(i in seq_along(A.species$values)) {
A.species$values[[i]][is.predominant] <- affs[[j]]$values[[i]][is.predominant]
}
Modified: pkg/CHNOSZ/R/util.plot.R
===================================================================
--- pkg/CHNOSZ/R/util.plot.R 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/R/util.plot.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -66,67 +66,66 @@
}
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") {
+ logaH2O=0, lty=2, lwd=1, col=par('fg'), xpoints=NULL, O2state="gas", plot.it=TRUE) {
# draw water stability limits
- # if we're on an Eh-pH diagram, or logfO2-pH diagram,
- # or logfO2-T or Eh-T
- # calculate them exactly (nicer looking lines), otherwise
- # (TODO) add them using affinity() and diagram()
-
- # get the x and y limits from the current plot
- pu <- par('usr')
- xlim <- pu[1:2]
- ylim <- pu[3:4]
- # exact lines
- # warning: Eh calculations are reliable only at a single T
+ # for Eh-pH, logfO2-pH, logfO2-T or Eh-T diagrams
+ # (i.e. redox variable is on the y axis)
+ 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('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)
- else if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
- else if(yaxis=="Eh") lines(xlim,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
- else if(yaxis=="pe") lines(xlim,convert(convert(logfO2,'E0',T=T,P=P,pH=xlim),"pe",T=T),lty=lty,lwd=lwd,col=col)
+ #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)
}
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)
- else if(yaxis=="Eh") lines(xlim,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
- else if(yaxis=="pe") lines(xlim,convert(convert(logfO2,'E0',T=T,P=P,pH=xlim),"pe",T=T),lty=lty,lwd=lwd,col=col)
+ #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)
}
} 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') {
- if(missing(T)) {
- xpoints <- seq(xlim[1],xlim[2],length.out=100)
- T <- envert(xpoints,"K")
- }
- }
- if(xaxis=='P') if(is.null(xpoints)) xpoints <- P
+ 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') lines(xpoints,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
- else lines(xpoints,logfO2,lty=lty,lwd=lwd,col=col)
+ 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') lines(xpoints,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
- else lines(xpoints,logfO2,lty=lty,lwd=lwd,col=col)
+ if(yaxis=='Eh') y.reduction <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
+ else y.reduction <- logfO2
}
- } else {
- # inexact lines
- #
}
+ if(yaxis=="Eh") y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=xpoints)
+ # 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)
+ }
+ # return the values
+ return(invisible(list(xpoints=xpoints, y.oxidation=y.oxidation, y.reduction=y.reduction)))
}
mtitle <- function(main, line=0, ...) {
Modified: pkg/CHNOSZ/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/demo/copper.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -38,7 +38,7 @@
# mosaic diagram with to speciate glycine as a function of pH
m <- mosaic(bases=Gly, pH=c(0, 16, 300), Eh=c(-0.6, 1.0, 300))
fill <- c(rep("lightgrey", 3), rep("white", 4), rep("lightblue", 4))
-d <- diagram(m$A.species, fill=fill, names=NULL, xaxs="i", yaxs="i")
+d <- diagram(m$A.species, fill=fill, names=NULL, xaxs="i", yaxs="i", fill.NA="pink2")
# to make the labels look nicer
names <- names[sort(unique(as.numeric(d$predominant)))]
for(i in 1:length(names)) {
@@ -54,7 +54,7 @@
}
# add glycine ionization lines
-d <- diagram(m$A.bases, add=TRUE, col="darkblue", lty=3, names=NULL)
+d <- diagram(m$A.bases, add=TRUE, col="darkblue", lty=3, names=NULL, limit.water=FALSE)
text(d$lx, -0.5, Gly, col="darkblue")
# add water lines and title
Modified: pkg/CHNOSZ/demo/mosaic.R
===================================================================
--- pkg/CHNOSZ/demo/mosaic.R 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/demo/mosaic.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -33,6 +33,7 @@
title(main=paste("Iron oxides, sulfides and carbonate in water, log(total S) = -6,",
"log(total C)=0, after Garrels and Christ, 1965", sep="\n"))
# overlay the carbonate basis species predominance fields
-diagram(m1$A.bases2, add=TRUE, col="blue", col.names="blue", lty=3)
+d <- diagram(m1$A.bases2, add=TRUE, col="blue", names=NULL, lty=3, limit.water=FALSE)
+text(d$lx, -0.8, as.expression(sapply(m1$A.bases2$species$name, expr.species)), col="blue")
# reset the database, as it was changed in this example
data(thermo)
Modified: pkg/CHNOSZ/demo/sources.R
===================================================================
--- pkg/CHNOSZ/demo/sources.R 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/demo/sources.R 2017-04-28 08:30:51 UTC (rev 190)
@@ -9,7 +9,7 @@
obigt.source <- unique(c(os1,os2))
obigt.source <- obigt.source[!is.na(obigt.source)]
# these all produce character(0) if the sources are all accounted for
-print("missing these sources (1) for thermodynamic properties:")
+print("missing these sources for thermodynamic properties:")
print(unique(obigt.source[!(obigt.source %in% ref.source)]))
# determine if all the reference sources are cited
# this should produce character(0)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/inst/NEWS 2017-04-28 08:30:51 UTC (rev 190)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-78 (2017-04-27)
+CHANGES IN CHNOSZ 1.0.8-79 (2017-04-28)
---------------------------------------
DOCUMENTATION:
@@ -58,6 +58,10 @@
- Add arguments `adj` and `dy` for x-alignment and y-offset of line
labels.
+- Add arguments `fill.NA` (color of empty areas) and `limit.water`
+ (assign NA to areas beyond water stability limits on Eh-pH and pe-pH
+ diagrams).
+
NEW FEATURES:
- Add ZC.col() for generating a red-grey-blue color scale from
@@ -161,7 +165,7 @@
the bug report and test.
- NaN values from equilibrate() are now preserved by diagram(),
- producing unlabeled fields rather than being mistakenly labeled with
+ producing empty (NA) fields rather than being mistakenly labeled with
the first species. Thanks to Grayson Boyer for the bug report.
OTHER CHANGES:
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/man/diagram.Rd 2017-04-28 08:30:51 UTC (rev 190)
@@ -14,7 +14,8 @@
ylog=TRUE, xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL,
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
lty=NULL, lwd=par("lwd"), dotted=NULL,
- col=par("col"), col.names=par("col"), fill=NULL, col.NA="black",
+ col=par("col"), col.names=par("col"), fill=NULL,
+ fill.NA="slategray1", limit.water=TRUE,
names=NULL, main=NULL, legend.x=NA, format.names=TRUE, adj=0.5, dy=0,
add=FALSE, plot.it=TRUE, tplot=TRUE, ...)
strip(affinity, ispecies = NULL, col = NULL, ns = NULL,
@@ -48,7 +49,8 @@
\item{col}{character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram), or colors of bars in a strip diagram (\code{strip})}
\item{col.names}{character, colors for labels of species}
\item{fill}{character, colors used to fill predominance fields}
- \item{col.NA}{character, color for grid points with NA values}
+ \item{fill.NA}{character, color for grid points with NA values}
+ \item{limit.water}{logical, set NA values beyond water stability limits?}
\item{names}{character, names of species for activity lines or predominance fields}
\item{main}{character, a main \code{\link{title}} for the plot; \code{NULL} means to plot no title}
\item{legend.x}{character, description of legend placement passed to \code{\link{legend}}}
@@ -89,6 +91,9 @@
\code{fill} determines the color of the predominance fields, \code{col} that of the boundary lines.
By default, \code{\link{heat.colors}} are used to fill the predominance fields in diagrams on the screen plot device.
\code{fill} can be any color specification, or the word \samp{rainbow}, \samp{heat}, \samp{terrain}, \samp{topo}, or \samp{cm}, indicating a palette from \pkg{grDevices}.
+\code{fill.NA} gives the color for empty fields, i.e. points for which NA values are present, possibly by using \code{\link{equilibrate}} at extreme conditions (see \code{test-diagram.Rd}).
+\code{fill.NA} is also used to specify the color outside the water stability limits on Eh-pH or pe-pH diagrams, when \code{limit.water} is TRUE.
+Note that the default for \code{fill.NA} is automatically changed to \samp{transparent} when \code{add} is TRUE.
As of CHNOSZ 1.0.8-11, a new default line-drawing procedure has been implemented.
This uses \code{\link{contour}} to draw smooth-looking diagonal and curved lines, at the expense of not coinciding exactly with the rectangular grid (which is still used for drawing colors).
@@ -210,9 +215,9 @@
"goethite", "melanterite", "pyrite"))
a <- affinity(pH=c(-1, 4, 256), pe=c(-5, 23, 256))
d <- diagram(a, main="Fe-S-O-H, after Majzlan et al., 2006")
-# the first four species show up along the top of the diagram
-stopifnot(all.equal(unique(t(d$predominant)[256,]), 1:4))
-water.lines(yaxis="pe")
+# 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)
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/util.plot.Rd
===================================================================
--- pkg/CHNOSZ/man/util.plot.Rd 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/man/util.plot.Rd 2017-04-28 08:30:51 UTC (rev 190)
@@ -27,7 +27,7 @@
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")
+ col = par("fg"), xpoints = NULL, O2state="gas", plot.it = TRUE)
mtitle(main, line=0, ...)
ZC.col(z)
}
@@ -64,6 +64,7 @@
\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}
\item{z}{numeric, set of values}
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-04-27 05:24:50 UTC (rev 189)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-04-28 08:30:51 UTC (rev 190)
@@ -244,7 +244,7 @@
as.chemical.formula(makeup(910))
```
-For organic species, a simple calculation of the average oxidation state of carbon (`r zc`) is possible given the species index, chemical formula, or elemental count:
+For organic species, a calculation of the average oxidation state of carbon (`r zc`) is possible given the species index, chemical formula, or elemental count:
```{r ZC_910, message=FALSE}
ZC(910)
ZC(info(910)$formula)
@@ -281,7 +281,7 @@
subcrt("water", T = c(400, 500, 600), P = c(200, 400, 600), grid = "P")$out$water
```
-```{r subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of H<sub>2</sub>O.", cache=TRUE, pngquant=pngquant, timeit=timeit}
+```{r subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of water.", cache=TRUE, pngquant=pngquant, timeit=timeit}
substuff <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
water <- substuff$out$water
plot(water$P, water$rho, type = "l")
@@ -545,8 +545,9 @@
```
```{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(-1, 1))
+a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
diagram(a, fill = "heat")
+water.lines()
```
Now we can calculate the affinities on an Eh-pH grid:
@@ -555,8 +556,10 @@
## Potential diagrams
-Given values of affinity, the <span style="color:green">`diagram()`</span> function uses the maximum affinity method to make a potential diagram (i.e. a Pourbaix diagram):
-```{r EhpH_plot, echo=2, eval=FALSE}
+Given values of affinity, the <span style="color:green">`diagram()`</span> function uses the maximum affinity method to make a potential diagram (i.e. a Pourbaix diagram).
+Areas corresponding to Eh-pH conditions beyond the stability limits of water are colored slate gray.
+A second function, <span style="color:green">`water.lines()`</span>, is used to draw lines at the water stability limits:
+```{r EhpH_plot, echo=-1, eval=FALSE}
```
Note that the calculation of affinity implies a non-equilibrium reference state of equal activities of species ([see above](#species-of-interest)).
@@ -572,7 +575,7 @@
The default colors for diagrams shown on the screen use R's `heat.colors()` palette.
Some arguments in <span style="color:green">`diagram()`</span> can be used to control the color, labels, and lines, and title.
The `tplot` argument turns off plot customizations used in CHNOSZ.
-Additional arguments are passed to R's plotting functions; here, we use `bty` to remove the box around the plot:
+Additional arguments are passed to R's plotting functions; here, we use `bty` to remove the box around the plot.
```{r EhpH_plot_color, echo=TRUE, eval=FALSE}
```
@@ -618,15 +621,16 @@
The key argument is `bases`, which identifies the candidate basis species, starting with the one in the current basis.
The other arguments, like those of <span style="color:green">`affinity()`</span>, specify the ranges of the variables; `res` indicates the grid resolution to use for each variable (the default is 128).
The first call to <span style="color:green">`diagram()`</span> plots the species of interest; the second adds the predominance fields of the basis species.
-Finally, <span style="color:green">`water.lines()`</span> is used to add the stability limits of water at the given temperature.
+We turn off the gray coloring beyond the water stability limits (`limit.water`) but plot the lines using <span style="color:green">`water.lines()`</span>:
-```{r copper_mosaic, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", message=FALSE, cache=TRUE, fig.cap="Copper minerals and aqueous complexes with chloride, 25 °C.", pngquant=pngquant, timeit=timeit}
+```{r copper_mosaic, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", message=FALSE, cache=TRUE, fig.cap="Copper minerals and aqueous complexes with chloride, 200 °C.", pngquant=pngquant, timeit=timeit}
T <- 200
res <- 200
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
m1 <- mosaic(bases, blend = TRUE, pH = c(0, 12, res), Eh=c(-1.2, 0.75, res), T=T)
-diagram(m1$A.species, lwd = 2, fill = NA)
-diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 2)
+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)
```
@@ -646,8 +650,10 @@
if (names(newvar) == "O2") basis("O2", "gas")
mosaicargs <- c(list(bases), blend=TRUE, pH=list(c(-2, 12, res)), newvar, T=T)
m1 <- do.call(mosaic, mosaicargs)
- diagram(m1$A.species, lwd = 2, fill = rev(topo.colors(10)))
- diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 3)
+ diagram(m1$A.species, lwd = 2, fill = rev(topo.colors(10)),
+ limit.water = FALSE)
+ diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 3,
+ limit.water = FALSE)
swap.basis(names(newvar), "e-")
}
par(mfrow = c(1, 3))
More information about the CHNOSZ-commits
mailing list