Refactor Pourbaix

Closes #2 (closed)

Will bump to version 0.3.11

R/Rphree_Pourbaix.R

 ### Licence: LGPL version 2.1
-## Time-stamp: "Last modified 2021-12-09 15:04:18 delucia"
+## Time-stamp: "Last modified 2022-07-01 17:07:49 delucia"
 
 ##' @title Revert a string
 ##' @param x character vector. All element of the vector will be
@@ -136,8 +136,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 plot logical, defaults to TRUE. Do you want to plot 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?
@@ -154,6 +152,10 @@ FormSelectedOutput <- function(sim, element) {
 ##'     Defaults to 1. Parallelization is achieved through DoParallel
 ##'     and foreach, and should work on both *nix and Windows using
 ##'     fork and sockets respectively
+##' @param offset numerical value, defaulting to 0, distance from the
+##'     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
 ##' @return invisibly returns a list containing all calculated
@@ -215,9 +217,9 @@ Pourbaix <- function(element,
                      ann.sub,
                      colors,
                      procs=1,
+                     offset=0,
                      plot=TRUE,
-                     ...)
-{  
+                     ...) {  
     phreeqc::phrLoadDatabaseString(db)
     phreeqc::phrSetOutputStringsOn(TRUE)
 
@@ -278,8 +280,8 @@ Pourbaix <- function(element,
 
     ## eliminate outside water stability region!
     ## TODO: check T/P dependance
-    ind_h2 <- which(combs$pe < - combs$pH)
-    ind_o2 <- which(combs$pe > 20.75 - combs$pH)
+    ind_h2 <- which(combs$pe < offset - combs$pH)
+    ind_o2 <- which(combs$pe > 20.75 - combs$pH - offset)
     ## we need to remember which simulation is suppressed to reform
     ## the grid afterwards
     unstable <- union(ind_o2,ind_h2)
@@ -329,7 +331,7 @@ Pourbaix <- function(element,
 
         ## distribute pe and pH into each chunk
         biginp <- lapply(ilist, function(x) {
-            inp <- Distribute(base2, "pe", x$pe);
+            inp <- Distribute(base2, "pe", x$pe)
             Distribute(inp, "pH", x$pH)
         })
         ## add "first" and "selout" blocks
@@ -470,7 +472,7 @@ Pourbaix <- function(element,
         ## 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){
+            if (!aqonly) {
                 ii <- which(is.aqueous==1)
                 ## we need also to check if ii is not empty!
                 if (length(ii)>0) {
@@ -484,7 +486,7 @@ 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)) {
@@ -505,7 +507,7 @@ Pourbaix <- function(element,
         F <- 96485.3365
         
         Eh <- seq(-1,1,0.25)
-        ppe <- Eh/2.3/8.3145/{Tc+273.15}*F
+        ppe <- Eh / 2.3 / 8.3145 / { Tc + 273.15 } * F
         ## Eh <-  2.3*8.3145*{Tc+273.15}* ppe/F
         ## par(new = TRUE)
         ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "")