[CHNOSZ-commits] r738 - in pkg/CHNOSZ: . R inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jul 23 16:39:52 CEST 2022
Author: jedick
Date: 2022-07-23 16:39:52 +0200 (Sat, 23 Jul 2022)
New Revision: 738
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/R/stack_mosaic.R
pkg/CHNOSZ/inst/NEWS.Rd
Log:
Add tests/test_stack_mosaic.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-07-22 23:52:24 UTC (rev 737)
+++ pkg/CHNOSZ/DESCRIPTION 2022-07-23 14:39:52 UTC (rev 738)
@@ -1,6 +1,6 @@
Date: 2022-07-23
Package: CHNOSZ
-Version: 1.9.9-29
+Version: 1.9.9-30
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2022-07-22 23:52:24 UTC (rev 737)
+++ pkg/CHNOSZ/R/solubility.R 2022-07-23 14:39:52 UTC (rev 738)
@@ -2,7 +2,7 @@
# 20181031 first version jmd
# 20181106 work on the output from affinity(); no "equilibrate()" needed!
# 20190117 add find.IS and test for dissociation reaction
-# 20210319 use list of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple minerals
+# 20210319 use vector of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple minerals
## if this file is interactively sourced, the following are also needed to provide unexported functions:
#source("equilibrate.R")
@@ -100,19 +100,19 @@
# The "nuts and bolts" of solubility calculations
# Moved from solubility() to solubility_calc() 20210318
solubility_calc <- function(aout, dissociate = NULL, find.IS = FALSE, in.terms.of = NULL) {
- ## concept: the logarithms of activities of species at equilibrium are equal to
+ ## Concept: the logarithms of activities of species at equilibrium are equal to
## Astar, the affinities calculated for unit activities of species
- ## is aout the output from mosaic() instead of affinity()?
+ ## Is aout the output from mosaic() instead of affinity()?
aout.save <- aout
thisfun <- aout$fun
if(thisfun=="mosaic") aout <- aout$A.species
- # get starting ionic strength (probably zero, but could be anything set by user)
+ # Get starting ionic strength (probably zero, but could be anything set by user)
IS <- aout$IS
if(is.null(IS)) IS <- 0
- ## to output loga.balance we need the balancing coefficients
+ ## To output loga.balance we need the balancing coefficients
bout <- balance(aout)
n.balance <- bout$n.balance
balance <- bout$balance
@@ -125,42 +125,42 @@
log10(a.balance)
}
- # for find.IS=TRUE, iterate to converge on ionic strength
+ # For find.IS = TRUE, iterate to converge on ionic strength
niter <- 1
while(TRUE) {
- ## the values in aout can be calculated for other than
+ ## The values in aout can be calculated for other than
## unit activities of species, so we have to take away the activites
Astar <- function(i) aout$values[[i]] + aout$species$logact[i]
loga.equil <- lapply(1:length(aout$values), Astar)
- ## for a dissociation on a *per reaction* (not system) basis, apply the divisor here
+ ## For a dissociation on a *per reaction* (not system) basis, apply the divisor here
## (can be used to reproduce Fig. 4 of Manning et al., 2013)
if(is.numeric(dissociate)) {
loga.equil <- lapply(loga.equil, "/", dissociate)
}
- # get the total activity of the balancing basis species
+ # Get the total activity of the balancing basis species
loga.balance <- logabfun(loga.equil, n.balance)
- # recalculate things for a dissociation reaction (like CaCO3 = Ca+2 + CO3+2)
+ # Recalculate things for a dissociation reaction (like CaCO3 = Ca+2 + CO3+2)
if(isTRUE(dissociate)) {
ndissoc <- 2
- # the number of dissociated products is the exponent in the activity product
+ # The number of dissociated products is the exponent in the activity product
loga.dissoc <- loga.balance / ndissoc
- # the contribution to affinity
+ # The contribution to affinity
Adissoc <- lapply(n.balance, "*", loga.dissoc)
- # adjust the affinity and get new equilibrium activities
+ # Adjust the affinity and get new equilibrium activities
aout$values <- mapply("-", aout$values, Adissoc, SIMPLIFY = FALSE)
loga.equil <- lapply(1:length(aout$values), Astar)
loga.balance <- logabfun(loga.equil, n.balance)
- # check that the new loga.balance == loga.dissoc
+ # Check that the new loga.balance == loga.dissoc
# TODO: does this work for non-1:1 species?
stopifnot(all.equal(loga.balance, loga.dissoc))
}
- # get the old IS
+ # Get the old IS
IS.old <- rep(IS, length.out = length(aout$values[[1]]))
- # calculate the ionic strength for unpaired ions: IS = sum(molality * Z^2) / 2
+ # Calculate the ionic strength for unpaired ions: IS = sum(molality * Z^2) / 2
sum.mZ2 <- 0
ion.names <- character()
# is there an ion in the basis species?
@@ -172,8 +172,8 @@
ion.names <- c(ion.names, basis.ion)
}
}
- # add ions present in the species of interest
- # instead of using aout$species$name, use info() to get formulas 20190309
+ # Add ions present in the species of interest
+ # Instead of using aout$species$name, use info() to get formulas 20190309
species.formulas <- suppressMessages(info(aout$species$ispecies, check.it = FALSE)$formula)
for(i in 1:length(loga.equil)) {
species.ion <- species.formulas[i]
@@ -184,15 +184,15 @@
}
}
IS <- sum.mZ2 / 2
- # report ions used in the ionic strength calculation
- if(find.IS & niter==1) message("solubility: ionic strength calculated for ", paste(ion.names, collapse=" "))
- # report the current ionic strength
- if(find.IS) message("solubility: (iteration ", niter, ") ionic strength range is ", paste(round(range(IS), 4), collapse=" "))
- # stop iterating if we reached the tolerance (or find.IS=FALSE)
+ # Report ions used in the ionic strength calculation
+ if(find.IS & niter == 1) message("solubility: ionic strength calculated for ", paste(ion.names, collapse = " "))
+ # Report the current ionic strength
+ if(find.IS) message("solubility: (iteration ", niter, ") ionic strength range is ", paste(round(range(IS), 4), collapse = " "))
+ # Stop iterating if we reached the tolerance (or find.IS = FALSE)
if(!find.IS | all(IS - IS.old < 1e-4)) break
- # on the first iteration, expand argument values for affinity() or mosaic()
+ # On the first iteration, expand argument values for affinity() or mosaic()
# (e.g. convert pH = c(0, 14) to pH = seq(0, 14, 256) so that it has the same length as the IS values)
- # we don't do this if aout$vals is NA 20190731
+ # We don't do this if aout$vals is NA 20190731
if(niter==1 & !all(is.na(aout$vals))) {
if(thisfun=="affinity") for(i in 1:length(aout$vals)) {
aout.save$args[[i]] <- aout$vals[[i]]
@@ -204,13 +204,13 @@
}
}
}
- # recalculate the affinity using the new IS
+ # Recalculate the affinity using the new IS
aout <- suppressMessages(do.call(thisfun, list(aout.save, IS = IS)))
if(thisfun=="mosaic") aout <- aout$A.species
niter <- niter + 1
}
- # do we want the solubility expressed in terms of
+ # Do we want the solubility expressed in terms of
# something other than the first basis species? 20190123
if(!is.null(in.terms.of)) {
# write the reaction between the basis species and the new species
@@ -221,12 +221,12 @@
loga.balance <- loga.balance - log10(coeff)
}
- # make the output (we don't deal with normalized formulas yet, so for now m.balance==n.balance)
- # indicate the function used to make this output 20190525
+ # Make the output (we don't deal with normalized formulas yet, so for now m.balance==n.balance)
+ # Indicate the function used to make this output 20190525
aout$fun <- "solubility"
- # add names to loga.equil 20190731
+ # Add names to loga.equil 20190731
names(loga.equil) <- aout$species$name
- c(aout, list(balance=bout$balance, m.balance=bout$n.balance, n.balance=bout$n.balance, in.terms.of=in.terms.of,
- loga.balance=loga.balance, Astar=loga.equil, loga.equil=loga.equil))
+ c(aout, list(balance = bout$balance, m.balance = bout$n.balance, n.balance = bout$n.balance, in.terms.of = in.terms.of,
+ loga.balance = loga.balance, Astar = loga.equil, loga.equil = loga.equil))
}
Modified: pkg/CHNOSZ/R/stack_mosaic.R
===================================================================
--- pkg/CHNOSZ/R/stack_mosaic.R 2022-07-22 23:52:24 UTC (rev 737)
+++ pkg/CHNOSZ/R/stack_mosaic.R 2022-07-23 14:39:52 UTC (rev 738)
@@ -17,7 +17,7 @@
# Set activities 20220722
if(!is.null(loga_aq)) {
iaq1 <- which(isp1$state == "aq")
- species(iaq1, loga_aq)
+ if(length(iaq1) > 0) species(iaq1, loga_aq)
}
# Calculate affinity of species1 while speciating bases (e.g. aqueous S species)
mosaic1 <- mosaic(bases, loga_aq = c(NA, loga_aq), ...)
@@ -30,11 +30,17 @@
# Set activities 20220722
if(!is.null(loga_aq)) {
iaq2 <- which(isp2$state == "aq")
- species(iaq2, loga_aq)
+ if(length(iaq2) > 0) species(iaq2, loga_aq)
}
+ # Use only predominant species1 as basis species (to speed up calculation) 20210224
+ predom <- diagram1$predominant
+ ipredom <- sort(unique(as.numeric(predom)))
+ for(i in seq_along(ipredom)) predom[diagram1$predominant == ipredom[i]] <- i
+ species1.predom <- species1[ipredom]
# Speciate bases again (NULL)
- # Take the predominant members of species1 (diagram1$predominant)
- mosaic2 <- mosaic(list(bases, species1), stable = list(NULL, diagram1$predominant), loga_aq = c(NA, loga_aq), ...)
+ # Take the predominant members of species1 (predom)
+ mosaic2 <- mosaic(list(bases, species1.predom), stable = list(NULL, predom), loga_aq = c(NA, loga_aq), ...)
+
# Set colors
col <- c(rep_len(col[[3]], length(species12)), rep_len(col[[2]], length(species2)))
col.names <- c(rep_len(col.names[[3]], length(species12)), rep_len(col.names[[2]], length(species2)))
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-07-22 23:52:24 UTC (rev 737)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-07-23 14:39:52 UTC (rev 738)
@@ -12,7 +12,7 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 1.9.9-29 (2022-07-23)}{
+\section{Changes in CHNOSZ version 1.9.9-30 (2022-07-23)}{
\subsection{MAJOR CHANGE}{
\itemize{
More information about the CHNOSZ-commits
mailing list