[CHNOSZ-commits] r939 - in pkg/CHNOSZ: . demo

[noreply at r-forge.r-project.org]

Author: jedick
Date: 2025-12-06 10:59:54 +0100 (Sat, 06 Dec 2025)
New Revision: 939

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/demo/NaCl.R
Log:
Update NaCl demo


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2025-11-11 12:56:06 UTC (rev 938)
+++ pkg/CHNOSZ/DESCRIPTION	2025-12-06 09:59:54 UTC (rev 939)
@@ -1,4 +1,4 @@
-Date: 2025-11-11
+Date: 2025-12-06
 Package: CHNOSZ
 Version: 2.2.0-6
 Title: Thermodynamic Calculations and Diagrams for Geochemistry

Modified: pkg/CHNOSZ/demo/NaCl.R
===================================================================
--- pkg/CHNOSZ/demo/NaCl.R	2025-11-11 12:56:06 UTC (rev 938)
+++ pkg/CHNOSZ/demo/NaCl.R	2025-12-06 09:59:54 UTC (rev 939)
@@ -10,21 +10,33 @@
 library(CHNOSZ)
 
 ## Uncomment these lines to make the plot with the g-function disabled
-#mod.OBIGT("Cl-", z=0)
-#mod.OBIGT("Na+", z=0)
+#mod.OBIGT("Cl-", z = 0)
+#mod.OBIGT("Na+", z = 0)
 
-# Start a new plot and show the experimental logK
-thermo.plot.new(xlim=c(0, 1000), ylim=c(-5.5, 1),
-  xlab=axis.label("T"), ylab=axis.label("logK"))
-expt <- read.csv(system.file("extdata/misc/SOJSH.csv", 
-  package="CHNOSZ"), as.is=TRUE)
-points(expt$T,expt$logK, pch=expt$pch)
+# Read the data and start a new plot
+expt <- read.csv(system.file("extdata/misc/SOJSH.csv", package = "CHNOSZ"), as.is = TRUE)
+expt$pch[expt$pch == 0] <- 15
+expt$pch[expt$pch == 1] <- 16
+expt$pch[expt$pch == 2] <- 17
+expt$pch[expt$pch == 5] <- 18
+thermo.plot.new(xlim = c(0, 1000), ylim = c(-5.5, 1), xlab = axis.label("T"), ylab = axis.label("logK"))
 
+# Use viridis palette
+expt_Ps <- unique(expt$P)
+cols <- hcl.colors(length(expt_Ps), palette = "Harmonic")
+# Loop over pressures
+for(i in 1:length(expt_Ps)) {
+  this_expt <- subset(expt, P == expt_Ps[i])
+  # Plot experimental logK
+  points(this_expt$T, this_expt$logK, cex = 1.5, pch = this_expt$pch, col = cols[i])
+}
+
 # We'll be at 9 distinct pressure conditions, including Psat
 # Psat is repeated to show "not considered" region
 # (T >= 355 degC; Fig. 6 of Shock et al., 1992)
-P <- c(list("Psat", "Psat"), as.list(seq(500, 4000, by=500)))
-# For each of those what's the range of temperature
+P <- c(list("Psat", "Psat"), as.list(seq(500, 4000, by = 500)))
+cols <- c(cols[1], cols)
+# The range of temperature for each pressure
 T <- list()
 T[[1]] <- seq(0, 354, 1)
 T[[2]] <- seq(354, 370, 1)
@@ -33,16 +45,17 @@
 T[[5]] <- seq(395, 920, 1)
 T[[6]] <- T[[7]] <- T[[8]] <- T[[9]] <- T[[10]] <- seq(400, 1000, 1)
 
-# Calculate and plot the logK
+# Define reaction for calculating logK
 species <- c("NaCl", "Na+", "Cl-")
 coeffs <- c(-1, 1, 1)
 logK <- numeric()
+# Loop over pressures
 for(i in 1:length(T)) {
-  s <- suppressWarnings(subcrt(species, coeffs, T=T[[i]], P=P[[i]]))
-  if(i==2) lty <- 3 else lty <- 1
-  lines(s$out$T, s$out$logK, lty=lty)
-  # keep the calculated values for each experimental condition (excluding Psat)
-  iexpt <- which(P[[i]]==expt$P)
+  s <- suppressWarnings(subcrt(species, coeffs, T = T[[i]], P = P[[i]]))
+  if(i == 2) lty <- 3 else lty <- 1
+  lines(s$out$T, s$out$logK, lwd = 2, lty = lty, col = cols[i])
+  # Keep the calculated values for each experimental condition (excluding Psat)
+  iexpt <- which(P[[i]] == expt$P)
   Texpt <- expt$T[iexpt]
   if(i > 2) logK <- c(logK, splinefun(s$out$T, s$out$logK)(Texpt))
 }
@@ -49,18 +62,17 @@
 
 # Add title, labels, and legends
 title(describe.reaction(s$reaction, states = 1))
-text(150, -0.1, quote(italic(P)[sat]), cex=1.2)
+text(150, -0.1, quote(italic(P)[sat]), cex = 1.2)
 text(462, -4, "500 bar")
 text(620, -4.3, "1000 bar")
 text(796, -4.3, "1500 bar")
-text(813, -1.4, "4000 bar")
-legend("bottomleft",pch=unique(expt$pch),
-  legend=c(unique(expt$source),tail(expt$source,1)), bty="n")
+text(813, -1.3, "4000 bar")
+legend("bottomleft", legend = unique(expt$source), pch = unique(expt$pch), pt.cex = c(1.2, 1.2, 1.4, 1.6), bty = "n")
 #mtitle(c(describe.reaction(s$reaction), expression(italic(P)[sat]~"and 500-4000 bar")))
 l1 <- quote("Revised HKF model with " * italic(g) * " function (Shock et al., 1992)")
 l2 <- "Non-recommended region (Shock et al., 1992, Fig. 6)"
-legend("topright", as.expression(c(l1, l2)), lty=c(1, 3), bty="n")
+legend("topright", as.expression(c(l1, l2)), lty = c(1, 3), bty = "n")
 
 # Test for average divergence (excluding Psat)
-expt <- expt[!expt$P %in% "Psat", ]
-stopifnot(mean(abs(logK - expt$logK)) < 0.09)
+expt_test <- expt[!expt$P %in% "Psat", ]
+stopifnot(mean(abs(logK - expt_test$logK)) < 0.09)