From d077a63509c3769b6ac34db8361517e1a6e1ab2e Mon Sep 17 00:00:00 2001 From: delucia Date: Wed, 13 Jul 2022 19:34:31 +0200 Subject: [PATCH 1/4] Bump to 0.3.11; syntax linter in Pourbaix() --- DESCRIPTION | 2 +- R/Rphree_Pourbaix.R | 21 ++++++++++++--------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7769949..ba49b60 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: RedModRphree Title: Utilities Leveraging the R Interface to the PHREEQC Geochemical Solver -Version: 0.3.10 +Version: 0.3.11 Authors@R: c(person(given = "Marco", family = "De Lucia", email = "delucia@gfz-potsdam.de", diff --git a/R/Rphree_Pourbaix.R b/R/Rphree_Pourbaix.R index c343937..dd2fb44 100644 --- a/R/Rphree_Pourbaix.R +++ b/R/Rphree_Pourbaix.R @@ -1,5 +1,6 @@ -### Licence: LGPL version 2.1 -## Time-stamp: "Last modified 2022-07-01 17:07:49 delucia" +## Functions to compute and plot Pourbaix diagrams + +## Time-stamp: "Last modified 2022-07-13 19:23:26 delucia" ##' @title Revert a string ##' @param x character vector. All element of the vector will be @@ -428,7 +429,9 @@ Pourbaix <- function(element, ## formulas of eqphases is.pp <- todisplay[is.aqueous==1] formulas[is.aqueous==1] <- res$SI$formula[match(is.pp, rownames(res$SI))] - } else eq <- NA + } else { + eq <- NA + } numvec <- match(vec, todisplay)## so that vec==todisplay[numvec] @@ -486,14 +489,14 @@ Pourbaix <- function(element, mymain <- paste("Pourbaix diagram for", ifelse(missing(element), "the specified system", ifelse(has_valence, valence, element)), "\n T=", Tc, "\u00B0C / P=", Patm, "atm") - } else + } else { mymain <- ann.title - + } if (missing(ann.sub)) { mysub <- paste("Molalities of dissolved elements:", paste(rownames(res$tot),res$tot[,1], collapse="; ")) - } else + } else { mysub <- ann.sub - + } title(main = mymain, sub = mysub) @@ -504,10 +507,10 @@ Pourbaix <- function(element, abline(0, -1, lty="dashed") abline(20.8, -1, lty="dashed") - F <- 96485.3365 + Faraday <- 96485.3365 Eh <- seq(-1,1,0.25) - ppe <- Eh / 2.3 / 8.3145 / { Tc + 273.15 } * F + ppe <- Faraday * Eh / 2.3 / 8.3145 / (Tc + 273.15) ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F ## par(new = TRUE) ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") -- GitLab From 2d24f7cf54a2314fb21bed9d3a053cc93a17b596 Mon Sep 17 00:00:00 2001 From: delucia Date: Wed, 13 Jul 2022 20:16:43 +0200 Subject: [PATCH 2/4] Tinkering with match.call, splitting plot --- R/Rphree_Pourbaix.R | 171 +++++++++++++++++++++++++------------------- 1 file changed, 97 insertions(+), 74 deletions(-) diff --git a/R/Rphree_Pourbaix.R b/R/Rphree_Pourbaix.R index dd2fb44..88f2280 100644 --- a/R/Rphree_Pourbaix.R +++ b/R/Rphree_Pourbaix.R @@ -1,6 +1,6 @@ ## Functions to compute and plot Pourbaix diagrams -## Time-stamp: "Last modified 2022-07-13 19:23:26 delucia" +## Time-stamp: "Last modified 2022-07-13 20:10:59 delucia" ##' @title Revert a string ##' @param x character vector. All element of the vector will be @@ -221,6 +221,10 @@ Pourbaix <- function(element, offset=0, plot=TRUE, ...) { + + ## store parameters + stored_pars <- as.list(match.call(expand.dots = FALSE)) + phreeqc::phrLoadDatabaseString(db) phreeqc::phrSetOutputStringsOn(TRUE) @@ -451,84 +455,103 @@ Pourbaix <- function(element, parsednames[is.aqueous==4] <- Raqlegexpr on.exit() - - ## plotting + ret <- list(mat=mat, + params=stored_pars, + displayed=todisplay, + formula=parsednames, + is.aq=is.aqueous, + phreeqc=list(input=c(first, selout), pqc=bigres, db=db, res=res), + dat=dat, + vec=vec, + aq=aq, + eq=eq, + levs=levs) + if (plot) { - if (missing(colors)) { - region_colors <- grDevices::terrain.colors(levs, alpha = 0.5) - names(region_colors) <- todisplay - } else if (!is.null(names(colors))) { - region_colors <- stats::na.omit(colors[match(todisplay, names(colors))]) - ## we now match the mat to the actual names of colors - numvec <- match(vec, names(region_colors))## so that vec==todisplay[numvec] - levs <- max(numvec, na.rm=TRUE) - mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) - } else { - region_colors <- colors - } - - - par(mar = c(5,5,5,5), las=1, family="serif") + PlotPourbaix(ret) + invisible(ret) + } + + return(ret) + +} + + + + +PlotPourbaix <- function(pbx, + aqonly=FALSE, + suppress, + ann.phases=TRUE, + ann.title, + ann.sub, + colors, + ...) { + + if (missing(colors)) { + region_colors <- grDevices::terrain.colors(levs, alpha = 0.5) + names(region_colors) <- todisplay + } else if (!is.null(names(colors))) { + region_colors <- stats::na.omit(colors[match(todisplay, names(colors))]) + ## we now match the mat to the actual names of colors + numvec <- match(vec, names(region_colors))## so that vec==todisplay[numvec] + levs <- max(numvec, na.rm=TRUE) + mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) + } else { + region_colors <- colors + } + + par(mar = c(5,5,5,5), las=1, family="serif") - image(x=pH, y=pe, z=mat, col = region_colors, useRaster=TRUE, axes=FALSE, ...) - - ## Do we want to plot the region names? - if (ann.phases) { - text(pH[centers[,1]], pe[centers[,2]], parsednames, font=is.aqueous, pos=3) - if (!aqonly) { - ii <- which(is.aqueous==1) - ## we need also to check if ii is not empty! - if (length(ii)>0) { - ## we plot these labels - text(pH[centers[ii,1]], pe[centers[ii,2]], Rformulas, font=3, pos=1) - } + image(x=pH, y=pe, z=mat, col = region_colors, useRaster=TRUE, axes=FALSE, ...) + + ## Do we want to plot the region names? + if (ann.phases) { + text(pH[centers[,1]], pe[centers[,2]], parsednames, font=is.aqueous, pos=3) + if (!aqonly) { + ii <- which(is.aqueous==1) + ## we need also to check if ii is not empty! + if (length(ii)>0) { + ## we plot these labels + text(pH[centers[ii,1]], pe[centers[ii,2]], Rformulas, font=3, pos=1) } } + } + + if (missing(ann.title)) { + mymain <- paste("Pourbaix diagram for", + ifelse(missing(element), "the specified system", ifelse(has_valence, valence, element)), + "\n T=", Tc, "\u00B0C / P=", Patm, "atm") + } else { + mymain <- ann.title + } + if (missing(ann.sub)) { + mysub <- paste("Molalities of dissolved elements:", paste(rownames(res$tot),res$tot[,1], collapse="; ")) + } else { + mysub <- ann.sub + } + + title(main = mymain, sub = mysub) + + axis(1, at=pretty(pH)) + axis(2, at=pretty(pe)) + + ## stability lines of water + abline(0, -1, lty="dashed") + abline(20.8, -1, lty="dashed") - if (missing(ann.title)) { - mymain <- paste("Pourbaix diagram for", - ifelse(missing(element), "the specified system", ifelse(has_valence, valence, element)), - "\n T=", Tc, "\u00B0C / P=", Patm, "atm") - } else { - mymain <- ann.title - } - if (missing(ann.sub)) { - mysub <- paste("Molalities of dissolved elements:", paste(rownames(res$tot),res$tot[,1], collapse="; ")) - } else { - mysub <- ann.sub - } - title(main = mymain, - sub = mysub) - - axis(1, at=pretty(pH)) - axis(2, at=pretty(pe)) - - ## stability lines of water - abline(0, -1, lty="dashed") - abline(20.8, -1, lty="dashed") - - Faraday <- 96485.3365 + Faraday <- 96485.3365 - Eh <- seq(-1,1,0.25) - ppe <- Faraday * Eh / 2.3 / 8.3145 / (Tc + 273.15) - ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F - ## par(new = TRUE) - ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") - axis(side=4, at = ppe, labels=Eh) - ## axis(3, at=1/{Tgrid+273.15}, labels=(Tgrid)) + Eh <- seq(-1,1,0.25) + ppe <- Faraday * Eh / 2.3 / 8.3145 / (Tc + 273.15) + ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F + ## par(new = TRUE) + ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") + axis(side=4, at = ppe, labels=Eh) + ## axis(3, at=1/{Tgrid+273.15}, labels=(Tgrid)) - mtext("Eh / V", side=4, line=3, las=0) - box() - } ## plotting - - invisible(list(mat=mat, - displayed=todisplay, - formula=parsednames, - is.aq=is.aqueous, - phreeqc=list(input=c(first, selout), pqc=bigres, db=db, res=res), - dat=dat, - vec=vec, - aq=aq, - eq=eq)) -} + mtext("Eh / V", side=4, line=3, las=0) + box() +} + -- GitLab From 649fcc4611788e324f557bf6f591a63716fbcde8 Mon Sep 17 00:00:00 2001 From: delucia Date: Thu, 14 Jul 2022 14:15:49 +0200 Subject: [PATCH 3/4] Splitted computation and plotting --- NAMESPACE | 1 + R/Rphree_Pourbaix.R | 180 ++++++++++++++++++++++++-------------------- man/PlotPourbaix.Rd | 35 +++++++++ man/Pourbaix.Rd | 42 +++-------- 4 files changed, 148 insertions(+), 110 deletions(-) create mode 100644 man/PlotPourbaix.Rd diff --git a/NAMESPACE b/NAMESPACE index be89e03..fbc9121 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,6 +36,7 @@ export(MatplotSingle) export(PQCRunAndDump) export(ParseFormula) export(PlotModsInSample) +export(PlotPourbaix) export(Pourbaix) export(RGetPhases) export(RPhreeExt) diff --git a/R/Rphree_Pourbaix.R b/R/Rphree_Pourbaix.R index 88f2280..5a0eb16 100644 --- a/R/Rphree_Pourbaix.R +++ b/R/Rphree_Pourbaix.R @@ -1,6 +1,6 @@ ## Functions to compute and plot Pourbaix diagrams -## Time-stamp: "Last modified 2022-07-13 20:10:59 delucia" +## Time-stamp: "Last modified 2022-07-14 14:07:34 delucia" ##' @title Revert a string ##' @param x character vector. All element of the vector will be @@ -110,14 +110,16 @@ FormSelectedOutput <- function(sim, element) { } -##' This function uses PHREEQC to calculate and plot a Pourbaix -##' diagram for the specified system at a given T and P. The -##' geochemical system is specified as a single PHREEQC input -##' specifying one "SOLUTION" while the actual values of pe and pH at -##' which the predominance shall be calculated must be specified by -##' the user at input. To try and prevent numerical instabilities, -##' which are quite frequent, the function only actually performs -##' PHREEQC calculations inside the water stability region. +##' This function uses PHREEQC to calculate (and optionally plot) a +##' Pourbaix diagram for the specified system at a given T and P. The +##' geochemical system is defined as a single PHREEQC input specifying +##' one "SOLUTION" while the actual values of pe and pH at which the +##' predominance shall be calculated must be specified by the user as +##' argument in the function call. To try and prevent numerical +##' instabilities, which are quite frequent, the function only +##' actually performs PHREEQC calculations inside the water stability +##' region. Through parameter \code{offset} the user can request a +##' given offset w.r.t. the limiting lines. ##' ##' @title Pourbaix (pe/pH) diagrams using PHREEQC ##' @param element character of length 1. The element for which the @@ -137,18 +139,6 @@ FormSelectedOutput <- function(sim, element) { ##' species are considered ##' @param suppress optional character vector. Species and minerals ##' contained here are not considered in the diagram -##' @param ann.phases Logical, defaults to TRUE. Do we write the -##' formulas of each phase in the (approximate) center of each -##' corresponding diagram region? -##' @param ann.title Optional string. If specified, it will replace -##' main title of the picture -##' @param ann.sub Optional string. If specified, it will replace -##' subtitle of the picture (usually this specifies the -##' concentrations) -##' @param colors Optional character vector containing the colors for -##' each region. If missing, grDevices::terrain.colors is used -##' internally for plotting. If vector is named, the names are -##' matched, if unnamed the colors are used sequentially ##' @param procs how many CPUs should be used for the computations? ##' Defaults to 1. Parallelization is achieved through DoParallel ##' and foreach, and should work on both *nix and Windows using @@ -157,8 +147,8 @@ FormSelectedOutput <- function(sim, element) { ##' limiting lines of water stability ##' @param plot logical, defaults to TRUE. Do you want to plot the ##' diagram? -##' @param ... further graphic parameters passed to the \code{image} -##' for plotting +##' @param ... further parameters passed to \code{PlotPourbaix} for +##' plotting ##' @return invisibly returns a list containing all calculated ##' simulations, input script for PHREEQC simulations, formulas, ##' names, etc @@ -213,17 +203,13 @@ Pourbaix <- function(element, db, aqonly=FALSE, suppress, - ann.phases=TRUE, - ann.title, - ann.sub, - colors, procs=1, offset=0, plot=TRUE, ...) { ## store parameters - stored_pars <- as.list(match.call(expand.dots = FALSE)) + store_call <- as.list(match.call(expand.dots = FALSE)) phreeqc::phrLoadDatabaseString(db) phreeqc::phrSetOutputStringsOn(TRUE) @@ -245,11 +231,14 @@ Pourbaix <- function(element, aa <- phreeqc::phrRunString(c(first, base)) con <- phreeqc::phrGetOutputStrings() + + has_element <- FALSE has_valence <- FALSE if (missing(element)) { utils::capture.output(res <- ReadOutVal(con)[[1]]) selout <- FormSelectedOutput(res) } else { + has_element <- TRUE ## first we distinguish if we want total amount of element or a specific valence if (length(grep("(", element, fixed=TRUE))>0) { ## now we call "valence" the full expression (including the parentheses) @@ -416,38 +405,38 @@ Pourbaix <- function(element, vec <- dat$Aqname } - todisplay <- stats::na.omit(unique(vec)) + displayed <- stats::na.omit(unique(vec)) if (!missing(suppress)) - todisplay <- todisplay[! todisplay %in% suppress] + displayed <- displayed[! displayed %in% suppress] - naqueous <- length(aq <- which(todisplay %in% unique(Aqmaxnam))) + naqueous <- length(aq <- which(displayed %in% unique(Aqmaxnam))) ## fonts 1 for eq phases, font 4 for aqueous - is.aqueous <- rep(1, length(todisplay)) - formulas <- rep("", length(todisplay)) + is.aqueous <- rep(1, length(displayed)) + formulas <- rep("", length(displayed)) is.aqueous[aq] <- 4 if (!aqonly) { - neqphases <- length(eq <- which(todisplay %in% unique(SImaxnam))) + neqphases <- length(eq <- which(displayed %in% unique(SImaxnam))) ## formulas of eqphases - is.pp <- todisplay[is.aqueous==1] + is.pp <- displayed[is.aqueous==1] formulas[is.aqueous==1] <- res$SI$formula[match(is.pp, rownames(res$SI))] } else { eq <- NA } - - numvec <- match(vec, todisplay)## so that vec==todisplay[numvec] - levs <- max(numvec, na.rm=TRUE) + molalities <- paste(rownames(res$tot), res$tot[,1], collapse="; ") + numvec <- match(vec, displayed)## so that vec==displayed[numvec] + levels <- max(numvec, na.rm=TRUE) mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) - centers <- matrix(0, ncol=2, nrow=length(todisplay)) - for (i in seq_along(todisplay)) { + centers <- matrix(0, ncol=2, nrow=length(displayed)) + for (i in seq_along(displayed)) { centers[i, ] <- floor(colMeans(which(mat == i, arr.ind = TRUE))) } ## Fix the names - parsednames <- todisplay + parsednames <- displayed legexpr <- FormulaToExpression(formulas[is.aqueous==1]) Rformulas <- parse(text = paste0("italic(", legexpr,")")) aqlegexpr <- FormulaToExpression(parsednames[is.aqueous==4]) @@ -455,103 +444,134 @@ Pourbaix <- function(element, parsednames[is.aqueous==4] <- Raqlegexpr on.exit() - ret <- list(mat=mat, - params=stored_pars, - displayed=todisplay, - formula=parsednames, - is.aq=is.aqueous, + + ret <- list(CALL=store_call, + pe=pe, + pH=pH, + Tc=Tc, + Patm=Patm, + mat=mat, + molalities=molalities, + centers=centers, + displayed=displayed, + parsednames=parsednames, + Rformulas=Rformulas, + is.aqueous=is.aqueous, + aqonly=aqonly, phreeqc=list(input=c(first, selout), pqc=bigres, db=db, res=res), - dat=dat, + #dat=dat, vec=vec, aq=aq, eq=eq, - levs=levs) + levels=levels, + suppr=suppr, + has_valence=has_valence, + has_element=has_element, + valence=ifelse(has_valence, valence, ""), + element=ifelse(has_element, element, "")) if (plot) { - PlotPourbaix(ret) + msg("Invoking PlotPourbaix") + PlotPourbaix(ret, ...) invisible(ret) + } else { + return(ret) } - - return(ret) - } - +##' @title Plot a Pourbaix diagram +##' @param pbx list returned by \code{Pourbaix()} +##' @param ann.phases Logical, defaults to TRUE. Do we write the +##' formulas of each phase in the (approximate) center of each +##' corresponding diagram region? +##' @param ann.title Optional string. If specified, it will replace +##' main title of the picture +##' @param ann.sub Optional string. If specified, it will replace +##' subtitle of the picture (usually this specifies the +##' concentrations) +##' @param colors Optional character vector containing the colors for +##' each region. If missing, grDevices::terrain.colors is used +##' internally for plotting. If vector is named, the names are +##' matched, if unnamed the colors are used sequentially +##' @param ... further parameters passed to \code{image} +##' @return +##' @author MDL +##' @export PlotPourbaix <- function(pbx, - aqonly=FALSE, - suppress, ann.phases=TRUE, ann.title, ann.sub, colors, ...) { + if (missing(colors)) { - region_colors <- grDevices::terrain.colors(levs, alpha = 0.5) - names(region_colors) <- todisplay + region_colors <- grDevices::terrain.colors(pbx$levels, alpha = 0.5) + names(region_colors) <- pbx$displayed } else if (!is.null(names(colors))) { - region_colors <- stats::na.omit(colors[match(todisplay, names(colors))]) + region_colors <- stats::na.omit(colors[match(pbx$displayed, names(colors))]) ## we now match the mat to the actual names of colors - numvec <- match(vec, names(region_colors))## so that vec==todisplay[numvec] - levs <- max(numvec, na.rm=TRUE) - mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) + numvec <- match(pbx$vec, names(region_colors))## so that vec==displayed[numvec] + mat <- matrix(numvec, ncol=length(pbx$pe), nrow=length(pbx$pH), byrow=TRUE) } else { region_colors <- colors } - + par(mar = c(5,5,5,5), las=1, family="serif") - - image(x=pH, y=pe, z=mat, col = region_colors, useRaster=TRUE, axes=FALSE, ...) - + + image(x=pbx$pH, y=pbx$pe, z=pbx$mat, col = region_colors, useRaster=TRUE, xlab="pH", ylab="pe", axes=FALSE, ...) + ## Do we want to plot the region names? if (ann.phases) { - text(pH[centers[,1]], pe[centers[,2]], parsednames, font=is.aqueous, pos=3) - if (!aqonly) { - ii <- which(is.aqueous==1) + text(pbx$pH[pbx$centers[,1]], pbx$pe[pbx$centers[,2]], pbx$parsednames, font=pbx$is.aqueous, pos=3) + if (!pbx$aqonly) { + ii <- which(pbx$is.aqueous==1) ## we need also to check if ii is not empty! if (length(ii)>0) { ## we plot these labels - text(pH[centers[ii,1]], pe[centers[ii,2]], Rformulas, font=3, pos=1) - } + text(pbx$pH[pbx$centers[ii,1]], pbx$pe[pbx$centers[ii,2]], pbx$Rformulas, font=3, pos=1) + } } } if (missing(ann.title)) { mymain <- paste("Pourbaix diagram for", - ifelse(missing(element), "the specified system", ifelse(has_valence, valence, element)), - "\n T=", Tc, "\u00B0C / P=", Patm, "atm") + ifelse(pbx$has_element, "the specified system", ifelse(pbx$has_valence, pbx$valence, pbx$element)), + "\n T=", pbx$Tc, "\u00B0C / P=", pbx$Patm, "atm") } else { mymain <- ann.title } if (missing(ann.sub)) { - mysub <- paste("Molalities of dissolved elements:", paste(rownames(res$tot),res$tot[,1], collapse="; ")) + mysub <- paste("Molalities of dissolved elements:", pbx$molalities) } else { mysub <- ann.sub } - + title(main = mymain, sub = mysub) - axis(1, at=pretty(pH)) - axis(2, at=pretty(pe)) + axis(1, at=pretty(pbx$pH)) + axis(2, at=pretty(pbx$pe)) ## stability lines of water abline(0, -1, lty="dashed") abline(20.8, -1, lty="dashed") - + Faraday <- 96485.3365 - + Eh <- seq(-1,1,0.25) - ppe <- Faraday * Eh / 2.3 / 8.3145 / (Tc + 273.15) + ppe <- Faraday * Eh / 2.3 / 8.3145 / (pbx$Tc + 273.15) ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F ## par(new = TRUE) ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") axis(side=4, at = ppe, labels=Eh) ## axis(3, at=1/{Tgrid+273.15}, labels=(Tgrid)) - + mtext("Eh / V", side=4, line=3, las=0) box() + + invisible() } diff --git a/man/PlotPourbaix.Rd b/man/PlotPourbaix.Rd new file mode 100644 index 0000000..4bb1621 --- /dev/null +++ b/man/PlotPourbaix.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Rphree_Pourbaix.R +\name{PlotPourbaix} +\alias{PlotPourbaix} +\title{Plot a Pourbaix diagram} +\usage{ +PlotPourbaix(pbx, ann.phases = TRUE, ann.title, ann.sub, colors, ...) +} +\arguments{ +\item{pbx}{list returned by \code{Pourbaix()}} + +\item{ann.phases}{Logical, defaults to TRUE. Do we write the +formulas of each phase in the (approximate) center of each +corresponding diagram region?} + +\item{ann.title}{Optional string. If specified, it will replace +main title of the picture} + +\item{ann.sub}{Optional string. If specified, it will replace +subtitle of the picture (usually this specifies the +concentrations)} + +\item{colors}{Optional character vector containing the colors for +each region. If missing, grDevices::terrain.colors is used +internally for plotting. If vector is named, the names are +matched, if unnamed the colors are used sequentially} + +\item{...}{further parameters passed to \code{image}} +} +\description{ +Plot a Pourbaix diagram +} +\author{ +MDL +} diff --git a/man/Pourbaix.Rd b/man/Pourbaix.Rd index 14d4ba1..5d4ba3e 100644 --- a/man/Pourbaix.Rd +++ b/man/Pourbaix.Rd @@ -15,10 +15,6 @@ Pourbaix( db, aqonly = FALSE, suppress, - ann.phases = TRUE, - ann.title, - ann.sub, - colors, procs = 1, offset = 0, plot = TRUE, @@ -53,22 +49,6 @@ species are considered} \item{suppress}{optional character vector. Species and minerals contained here are not considered in the diagram} -\item{ann.phases}{Logical, defaults to TRUE. Do we write the -formulas of each phase in the (approximate) center of each -corresponding diagram region?} - -\item{ann.title}{Optional string. If specified, it will replace -main title of the picture} - -\item{ann.sub}{Optional string. If specified, it will replace -subtitle of the picture (usually this specifies the -concentrations)} - -\item{colors}{Optional character vector containing the colors for -each region. If missing, grDevices::terrain.colors is used -internally for plotting. If vector is named, the names are -matched, if unnamed the colors are used sequentially} - \item{procs}{how many CPUs should be used for the computations? Defaults to 1. Parallelization is achieved through DoParallel and foreach, and should work on both *nix and Windows using @@ -80,8 +60,8 @@ limiting lines of water stability} \item{plot}{logical, defaults to TRUE. Do you want to plot the diagram?} -\item{...}{further graphic parameters passed to the \code{image} -for plotting} +\item{...}{further parameters passed to \code{PlotPourbaix} for +plotting} } \value{ invisibly returns a list containing all calculated @@ -89,14 +69,16 @@ invisibly returns a list containing all calculated names, etc } \description{ -This function uses PHREEQC to calculate and plot a Pourbaix -diagram for the specified system at a given T and P. The -geochemical system is specified as a single PHREEQC input -specifying one "SOLUTION" while the actual values of pe and pH at -which the predominance shall be calculated must be specified by -the user at input. To try and prevent numerical instabilities, -which are quite frequent, the function only actually performs -PHREEQC calculations inside the water stability region. +This function uses PHREEQC to calculate (and optionally plot) a +Pourbaix diagram for the specified system at a given T and P. The +geochemical system is defined as a single PHREEQC input specifying +one "SOLUTION" while the actual values of pe and pH at which the +predominance shall be calculated must be specified by the user as +argument in the function call. To try and prevent numerical +instabilities, which are quite frequent, the function only +actually performs PHREEQC calculations inside the water stability +region. Through parameter \code{offset} the user can request a +given offset w.r.t. the limiting lines. } \examples{ base <- c("SOLUTION 1", -- GitLab From 3ffef1367a67d5968a49afc0c13b7982af70e9df Mon Sep 17 00:00:00 2001 From: delucia Date: Thu, 14 Jul 2022 15:56:27 +0200 Subject: [PATCH 4/4] Some docu improvements in PlotPourbaix --- R/Rphree_Pourbaix.R | 39 +++++++++++++++++++++++++++++++-------- man/PlotPourbaix.Rd | 28 ++++++++++++++++++++++++++++ man/Pourbaix.Rd | 5 +++-- 3 files changed, 62 insertions(+), 10 deletions(-) diff --git a/R/Rphree_Pourbaix.R b/R/Rphree_Pourbaix.R index 5a0eb16..9b5f968 100644 --- a/R/Rphree_Pourbaix.R +++ b/R/Rphree_Pourbaix.R @@ -1,6 +1,6 @@ ## Functions to compute and plot Pourbaix diagrams -## Time-stamp: "Last modified 2022-07-14 14:07:34 delucia" +## Time-stamp: "Last modified 2022-07-14 15:55:47 delucia" ##' @title Revert a string ##' @param x character vector. All element of the vector will be @@ -106,7 +106,7 @@ FormSelectedOutput <- function(sim, element) { " -percent_error true", Totals, Species, SI) - return(SelOut) + SelOut } @@ -171,8 +171,9 @@ FormSelectedOutput <- function(sim, element) { ##' par(mfrow=c(1,2)) ##' a <- Pourbaix(element="Ca", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, ##' Patm=1, base=base, first="", ann.title="Custom title", db=llnl.dat) -##' b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, -##' Patm=1, base=base, first="", ann.sub="Custom subtitle", db=llnl.dat) +##' b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), +##' Tc= 25, Patm=1, base=base, first="", ann.sub="Custom subtitle", +##' db=llnl.dat, offset=0.2) ##' ##' ## Now we restrict the speciation to a fixed valence state of the element, ##' ## here Fe(2) and Cu(2), only aqueous species: @@ -475,7 +476,7 @@ Pourbaix <- function(element, PlotPourbaix(ret, ...) invisible(ret) } else { - return(ret) + ret } } @@ -496,17 +497,39 @@ Pourbaix <- function(element, ##' internally for plotting. If vector is named, the names are ##' matched, if unnamed the colors are used sequentially ##' @param ... further parameters passed to \code{image} -##' @return +##' @return nothing ##' @author MDL ##' @export +##' @examples +##' base <- c("SOLUTION 1", +##' "units mol/kgw", +##' "temp 25.0", +##' "pressure 1", +##' "-water 1", +##' "pH 7", +##' "pe 4.0", +##' "Ca 0.025", +##' "Cu 0.001", +##' "Cl 0.05122 charge", +##' "Fe 3.0E-04", +##' "END") +##' +##' b <- Pourbaix(pe=seq(-10,12, length=51), pH=seq(2,12, length=51), +##' Patm=1, Tc=25, base=base, db=llnl.dat, plot=FALSE) +##' +##' mycol <- c("Cu"="grey", "Cl-"="white","Hematite"="light blue", "Nantokite"="red", +##' "Delafossite"="orange", "Ferrite-Cu"="gold", "Magnetite"="light green", +##' "Atacamite"="grey90") +##' +##' par(mfrow=c(1,2)) +##' PlotPourbaix(b, ann.title = "Default colors") +##' PlotPourbaix(b, colors = mycol, ann.title = "Custom colors via named vector") PlotPourbaix <- function(pbx, ann.phases=TRUE, ann.title, ann.sub, colors, ...) { - - if (missing(colors)) { region_colors <- grDevices::terrain.colors(pbx$levels, alpha = 0.5) names(region_colors) <- pbx$displayed diff --git a/man/PlotPourbaix.Rd b/man/PlotPourbaix.Rd index 4bb1621..5d4f50e 100644 --- a/man/PlotPourbaix.Rd +++ b/man/PlotPourbaix.Rd @@ -27,9 +27,37 @@ matched, if unnamed the colors are used sequentially} \item{...}{further parameters passed to \code{image}} } +\value{ +nothing +} \description{ Plot a Pourbaix diagram } +\examples{ +base <- c("SOLUTION 1", + "units mol/kgw", + "temp 25.0", + "pressure 1", + "-water 1", + "pH 7", + "pe 4.0", + "Ca 0.025", + "Cu 0.001", + "Cl 0.05122 charge", + "Fe 3.0E-04", + "END") + +b <- Pourbaix(pe=seq(-10,12, length=51), pH=seq(2,12, length=51), + Patm=1, Tc=25, base=base, db=llnl.dat, plot=FALSE) + +mycol <- c("Cu"="grey", "Cl-"="white","Hematite"="light blue", "Nantokite"="red", + "Delafossite"="orange", "Ferrite-Cu"="gold", "Magnetite"="light green", + "Atacamite"="grey90") + +par(mfrow=c(1,2)) +PlotPourbaix(b, ann.title = "Default colors") +PlotPourbaix(b, colors = mycol, ann.title = "Custom colors via named vector") +} \author{ MDL } diff --git a/man/Pourbaix.Rd b/man/Pourbaix.Rd index 5d4ba3e..7c0219c 100644 --- a/man/Pourbaix.Rd +++ b/man/Pourbaix.Rd @@ -97,8 +97,9 @@ base <- c("SOLUTION 1", par(mfrow=c(1,2)) a <- Pourbaix(element="Ca", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, Patm=1, base=base, first="", ann.title="Custom title", db=llnl.dat) -b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, - Patm=1, base=base, first="", ann.sub="Custom subtitle", db=llnl.dat) +b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), + Tc= 25, Patm=1, base=base, first="", ann.sub="Custom subtitle", + db=llnl.dat, offset=0.2) ## Now we restrict the speciation to a fixed valence state of the element, ## here Fe(2) and Cu(2), only aqueous species: -- GitLab