Nothing
inst/doc/eos-regress.R
In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----digits6, include=FALSE---------------------------------------------------
options(digits = 6)
## ----HTML, include=FALSE------------------------------------------------------
## Some frequently used HTML expressions
V0 <- "<i>V</i>°"
Cp0 <- "<i>C<sub>P</sub></i>°"
c1 <- "<i>c</i><sub>1</sub>"
c2 <- "<i>c</i><sub>2</sub>"
a1 <- "<i>a</i><sub>1</sub>"
a2 <- "<i>a</i><sub>2</sub>"
a3 <- "<i>a</i><sub>3</sub>"
a4 <- "<i>a</i><sub>4</sub>"
h4sio4 <- "H<sub>4</sub>SiO<sub>4</sub>"
sio2 <- "SiO<sub>2</sub>"
h2o <- "H<sub>2</sub>O"
ch4 <- "CH<sub>4</sub>"
wPrTr <- "ω<sub><i>P<sub>r</sub></i>,<i>T<sub>r</sub></i></sub>"
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
# Invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
# Adjust plot margins
knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .3, .3)) # smaller margin on top and right
})
# Use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) {
pngquant <- NULL
} else {
pngquant <- "--speed=1 --quality=0-50"
}
## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
)
## ----library_CHNOSZ-----------------------------------------------------------
library(CHNOSZ)
reset()
## ----Cpdat--------------------------------------------------------------------
file <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
Cpdat <- read.csv(file)
# Use data for CH4
Cpdat <- Cpdat[Cpdat$species == "CH4", ]
Cpdat <- Cpdat[, -1]
Cpdat$P <- convert(Cpdat$P, "bar")
## ----EOSregress---------------------------------------------------------------
var <- c("invTTheta2", "TXBorn")
Cplm_high <- EOSregress(Cpdat, var)
Cplm_low <- EOSregress(Cpdat, var, T.max = 600)
Cplm_low$coefficients
## ----EOSplot, fig.margin = TRUE, fig.cap = "Heat capacity of aqueous methane.", fig.width=3.5, fig.height=3.5, cache=TRUE, results="hide", message=FALSE, echo=FALSE, out.width=672, out.height=336, pngquant=pngquant----
EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))
EOSplot(Cpdat, coeficients = Cplm_high, add = TRUE, lty = 3)
PS01_data <- convert(EOScoeffs("CH4", "Cp"), "J")
EOSplot(Cpdat, coefficients = PS01_data, add = TRUE, lty = 2, col = "blue")
## ----EOSplot, eval=FALSE------------------------------------------------------
# EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))
# EOSplot(Cpdat, coeficients = Cplm_high, add = TRUE, lty = 3)
# PS01_data <- convert(EOScoeffs("CH4", "Cp"), "J")
# EOSplot(Cpdat, coefficients = PS01_data, add = TRUE, lty = 2, col = "blue")
## ----Cpcoeffs, message=FALSE--------------------------------------------------
convert(EOScoeffs("CH4", "Cp"), "J")
## ----Vdat---------------------------------------------------------------------
file <- system.file("extdata/cpetc/HWM96_V.csv", package = "CHNOSZ")
Vdat <- read.csv(file)
# Use data for CH4 near 280 bar
Vdat <- Vdat[Vdat$species == "CH4", ]
Vdat <- Vdat[abs(Vdat$P - 28) < 0.1, ]
Vdat <- Vdat[, -1]
Vdat$P <- convert(Vdat$P, "bar")
## ----Vdat_non-----------------------------------------------------------------
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm_low$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat_non <- data.frame(T = Vdat$T, P = Vdat$P, V = V_non)
## ----Vdat_non_regress, message=FALSE------------------------------------------
var <- "invTTheta"
Vnonlm <- EOSregress(Vdat_non, var, T.max = 450)
Vcoeffs <- round(c(Vnonlm$coefficients, QBorn = Vomega), 1)
Vcoeffs_database <- convert(EOScoeffs("CH4", "V"), "J")
## ----Vplot, fig.margin=TRUE, results="hide", message=FALSE, echo=FALSE, fig.width=3.5, fig.height=7, fig.cap="Volume of aqueous methane.", out.width=672, out.height=672, pngquant=pngquant----
par(mfrow = c(2, 1))
# plot 1
EOSplot(Vdat, coefficients = Vcoeffs)
EOSplot(Vdat, coefficients = Vcoeffs_database, add = TRUE, lty = 2)
# plot 2
EOSplot(Vdat, coefficients = Vcoeffs_database, T.plot = 600, lty = 2)
EOSplot(Vdat, coefficients = Vcoeffs, add = TRUE)
## ----Vplot, eval=FALSE--------------------------------------------------------
# par(mfrow = c(2, 1))
# # plot 1
# EOSplot(Vdat, coefficients = Vcoeffs)
# EOSplot(Vdat, coefficients = Vcoeffs_database, add = TRUE, lty = 2)
# # plot 2
# EOSplot(Vdat, coefficients = Vcoeffs_database, T.plot = 600, lty = 2)
# EOSplot(Vdat, coefficients = Vcoeffs, add = TRUE)
## ----Nadat--------------------------------------------------------------------
T <- convert(seq(0, 600, 50), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
Nadat <- prop.PT[, c("T", "P", "Cp")]
## ----Nalm, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na<sup>+</sup> (inapplicable: constant ω).", out.width=672, out.height=336, pngquant=pngquant----
var <- c("invTTheta2", "TXBorn")
Nalm <- EOSregress(Nadat, var, T.max = 600)
EOSplot(Nadat, coefficients = Nalm$coefficients, fun.legend = NULL)
EOSplot(Nadat, add = TRUE, lty = 3)
## ----Navars1------------------------------------------------------------------
var1 <- c("invTTheta2", "Cp_s_var")
omega.guess <- coef(Nalm)[3]
## ----Nawhile, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na<sup>+</sup> (variable ω).", out.width=672, out.height=336, pngquant=pngquant----
diff.omega <- 999
while(abs(diff.omega) > 1) {
Nalm1 <- EOSregress(Nadat, var1, omega.PrTr = tail(omega.guess, 1), Z = 1)
omega.guess <- c(omega.guess, coef(Nalm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(Nadat, coefficients = signif(coef(Nalm1), 6),
omega.PrTr = tail(omega.guess, 1), Z = 1)
convert(EOScoeffs("Na+", "Cp"), "J")
## ----NaVolume, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Volume of Na<sup>+</sup> (variable ω).", results="hide", message=FALSE, echo=FALSE, out.width=672, out.height=336, pngquant=pngquant----
T <- convert(seq(0, 600, 25), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
NaVdat <- prop.PT[, c("T", "P", "V")]
var1 <- c("invTTheta", "V_s_var")
omega.guess <- 1400000
diff.omega <- 999
while(abs(diff.omega) > 1) {
NaVlm1 <- EOSregress(NaVdat, var1,
omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1)
omega.guess <- c(omega.guess, coef(NaVlm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(NaVdat, coefficients = signif(coef(NaVlm1), 6),
omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1,
fun.legend = "bottomleft")
coefficients <- convert(EOScoeffs("Na+", "V", P = 1000), "J")
names(coefficients)[3] <- "V_s_var"
EOSplot(NaVdat, coefficients = coefficients, Z = 1, add = TRUE, lty = 2,
omega.PrTr = convert(coefficients["V_s_var"], "joules"))
## ----NaVolume, eval=FALSE-----------------------------------------------------
# T <- convert(seq(0, 600, 25), "K")
# P <- 1000
# prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
# NaVdat <- prop.PT[, c("T", "P", "V")]
# var1 <- c("invTTheta", "V_s_var")
# omega.guess <- 1400000
# diff.omega <- 999
# while(abs(diff.omega) > 1) {
# NaVlm1 <- EOSregress(NaVdat, var1,
# omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1)
# omega.guess <- c(omega.guess, coef(NaVlm1)[3])
# diff.omega <- tail(diff(omega.guess), 1)
# }
# EOSplot(NaVdat, coefficients = signif(coef(NaVlm1), 6),
# omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1,
# fun.legend = "bottomleft")
# coefficients <- convert(EOScoeffs("Na+", "V", P = 1000), "J")
# names(coefficients)[3] <- "V_s_var"
# EOSplot(NaVdat, coefficients = coefficients, Z = 1, add = TRUE, lty = 2,
# omega.PrTr = convert(coefficients["V_s_var"], "joules"))
## ----AS04, message=FALSE------------------------------------------------------
add.OBIGT("AS04")
## ----SiO2_2H2O, message=FALSE-------------------------------------------------
s_25C <- subcrt(c("SiO2", "H2O"), c(1, 2), T = 25)$out
s_near25 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(20, 30, length.out=50))$out
s_lowT <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 100, 10))$out
s_Psat <- subcrt(c("SiO2", "H2O"), c(1, 2))$out
s_P500 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 500)$out
s_P1000 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 1000)$out
## ----new_H4SiO4---------------------------------------------------------------
mod.OBIGT("calc-H4SiO4", formula = "H4SiO4", ref1 = "this_vignette",
date = as.character(Sys.Date()), G = s_25C$G, H = s_25C$H, S = s_25C$S,
Cp = s_25C$Cp, V = s_25C$V, z = 0)
## ----substuff-----------------------------------------------------------------
substuff <- rbind(s_near25, s_lowT, s_Psat, s_P500, s_P1000)
substuff$T <- convert(substuff$T, "K")
## ----Cp_H4SiO4, results="hide"------------------------------------------------
Cpdat <- substuff[, c("T", "P", "Cp")]
var <- c("invTTheta2", "TXBorn")
Cplm <- EOSregress(Cpdat, var)
Cpcoeffs <- Cplm$coefficients
mod.OBIGT("calc-H4SiO4", c1 = Cpcoeffs[1],
c2 = Cpcoeffs[2]/10000, omega = Cpcoeffs[3]/100000)
## ----V_H4SiO4_nonsolvation----------------------------------------------------
Vdat <- substuff[, c("T", "P", "V")]
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat$V <- V_non
## ----V_H4SiO4, results="hide"-------------------------------------------------
var <- c("invPPsi", "invTTheta", "invPPsiTTheta")
Vlm <- EOSregress(Vdat, var)
Vcoeffs <- convert(Vlm$coefficients, "joules")
mod.OBIGT("calc-H4SiO4", a1 = Vcoeffs[1]*10, a2 = Vcoeffs[2]/100,
a3 = Vcoeffs[3], a4 = Vcoeffs[4]/10000)
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----info_H4SiO4, message=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------
info(info(c("calc-H4SiO4", "H4SiO4")))
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----subcrt_H4SiO4, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=FALSE, results="hide", message=FALSE, out.width="100%", cache=TRUE, fig.cap="Comparison of H<sub>4</sub>SiO<sub>4</sub> pseudospecies.", pngquant=pngquant----
s1 <- subcrt(c("calc-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
plot(s1$out$T, s1$out$G, type = "l", ylim = c(-500, 2000),
xlab = axis.label("T"), ylab = axis.label("DG0"))
s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
lines(s2$out$T, s2$out$G, lty = 2)
abline(h = 0, lty = 3)
legend("topright", legend = c("calc-H4SiO4 (this vignette)",
"H4SiO4 (Stef\u00e1nsson, 2001)"), lty = c(1, 2), bty = "n")
text(225, 1000, describe.reaction(s1$reaction))
## ----subcrt_H4SiO4, eval=FALSE------------------------------------------------
# s1 <- subcrt(c("calc-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
# plot(s1$out$T, s1$out$G, type = "l", ylim = c(-500, 2000),
# xlab = axis.label("T"), ylab = axis.label("DG0"))
# s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
# lines(s2$out$T, s2$out$G, lty = 2)
# abline(h = 0, lty = 3)
# legend("topright", legend = c("calc-H4SiO4 (this vignette)",
# "H4SiO4 (Stef\u00e1nsson, 2001)"), lty = c(1, 2), bty = "n")
# text(225, 1000, describe.reaction(s1$reaction))
## ----activity_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=TRUE, results="hide", message=FALSE, out.width="100%", cache=TRUE, fig.cap="Activity diagram for K<sub>2</sub>O-Al<sub>2</sub>O<sub>3</sub>-SiO<sub>2</sub>-H<sub>2</sub>O.", pngquant=pngquant----
basis(c("Al+3", "H4SiO4", "K+", "H2O", "H+", "O2"))
species(c("gibbsite", "muscovite", "kaolinite", "pyrophyllite", "K-feldspar"))
a <- affinity(H4SiO4 = c(-8, 0, 300), `K+` = c(-1, 8, 300))
diagram(a, ylab = ratlab("K+"), fill = "terrain", yline = 1.7)
legend("bottomleft", describe.property(c("T", "P"), c(25, 1)), bty = "n")
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
CHNOSZ documentation built on May 29, 2024, 3:30 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----digits6, include=FALSE---------------------------------------------------
options(digits = 6)
## ----HTML, include=FALSE------------------------------------------------------
## Some frequently used HTML expressions
V0 <- "<i>V</i>°"
Cp0 <- "<i>C<sub>P</sub></i>°"
c1 <- "<i>c</i><sub>1</sub>"
c2 <- "<i>c</i><sub>2</sub>"
a1 <- "<i>a</i><sub>1</sub>"
a2 <- "<i>a</i><sub>2</sub>"
a3 <- "<i>a</i><sub>3</sub>"
a4 <- "<i>a</i><sub>4</sub>"
h4sio4 <- "H<sub>4</sub>SiO<sub>4</sub>"
sio2 <- "SiO<sub>2</sub>"
h2o <- "H<sub>2</sub>O"
ch4 <- "CH<sub>4</sub>"
wPrTr <- "ω<sub><i>P<sub>r</sub></i>,<i>T<sub>r</sub></i></sub>"
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
# Invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
# Adjust plot margins
knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .3, .3)) # smaller margin on top and right
})
# Use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) {
pngquant <- NULL
} else {
pngquant <- "--speed=1 --quality=0-50"
}
## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
)
## ----library_CHNOSZ-----------------------------------------------------------
library(CHNOSZ)
reset()
## ----Cpdat--------------------------------------------------------------------
file <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
Cpdat <- read.csv(file)
# Use data for CH4
Cpdat <- Cpdat[Cpdat$species == "CH4", ]
Cpdat <- Cpdat[, -1]
Cpdat$P <- convert(Cpdat$P, "bar")
## ----EOSregress---------------------------------------------------------------
var <- c("invTTheta2", "TXBorn")
Cplm_high <- EOSregress(Cpdat, var)
Cplm_low <- EOSregress(Cpdat, var, T.max = 600)
Cplm_low$coefficients
## ----EOSplot, fig.margin = TRUE, fig.cap = "Heat capacity of aqueous methane.", fig.width=3.5, fig.height=3.5, cache=TRUE, results="hide", message=FALSE, echo=FALSE, out.width=672, out.height=336, pngquant=pngquant----
EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))
EOSplot(Cpdat, coeficients = Cplm_high, add = TRUE, lty = 3)
PS01_data <- convert(EOScoeffs("CH4", "Cp"), "J")
EOSplot(Cpdat, coefficients = PS01_data, add = TRUE, lty = 2, col = "blue")
## ----EOSplot, eval=FALSE------------------------------------------------------
# EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))
# EOSplot(Cpdat, coeficients = Cplm_high, add = TRUE, lty = 3)
# PS01_data <- convert(EOScoeffs("CH4", "Cp"), "J")
# EOSplot(Cpdat, coefficients = PS01_data, add = TRUE, lty = 2, col = "blue")
## ----Cpcoeffs, message=FALSE--------------------------------------------------
convert(EOScoeffs("CH4", "Cp"), "J")
## ----Vdat---------------------------------------------------------------------
file <- system.file("extdata/cpetc/HWM96_V.csv", package = "CHNOSZ")
Vdat <- read.csv(file)
# Use data for CH4 near 280 bar
Vdat <- Vdat[Vdat$species == "CH4", ]
Vdat <- Vdat[abs(Vdat$P - 28) < 0.1, ]
Vdat <- Vdat[, -1]
Vdat$P <- convert(Vdat$P, "bar")
## ----Vdat_non-----------------------------------------------------------------
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm_low$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat_non <- data.frame(T = Vdat$T, P = Vdat$P, V = V_non)
## ----Vdat_non_regress, message=FALSE------------------------------------------
var <- "invTTheta"
Vnonlm <- EOSregress(Vdat_non, var, T.max = 450)
Vcoeffs <- round(c(Vnonlm$coefficients, QBorn = Vomega), 1)
Vcoeffs_database <- convert(EOScoeffs("CH4", "V"), "J")
## ----Vplot, fig.margin=TRUE, results="hide", message=FALSE, echo=FALSE, fig.width=3.5, fig.height=7, fig.cap="Volume of aqueous methane.", out.width=672, out.height=672, pngquant=pngquant----
par(mfrow = c(2, 1))
# plot 1
EOSplot(Vdat, coefficients = Vcoeffs)
EOSplot(Vdat, coefficients = Vcoeffs_database, add = TRUE, lty = 2)
# plot 2
EOSplot(Vdat, coefficients = Vcoeffs_database, T.plot = 600, lty = 2)
EOSplot(Vdat, coefficients = Vcoeffs, add = TRUE)
## ----Vplot, eval=FALSE--------------------------------------------------------
# par(mfrow = c(2, 1))
# # plot 1
# EOSplot(Vdat, coefficients = Vcoeffs)
# EOSplot(Vdat, coefficients = Vcoeffs_database, add = TRUE, lty = 2)
# # plot 2
# EOSplot(Vdat, coefficients = Vcoeffs_database, T.plot = 600, lty = 2)
# EOSplot(Vdat, coefficients = Vcoeffs, add = TRUE)
## ----Nadat--------------------------------------------------------------------
T <- convert(seq(0, 600, 50), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
Nadat <- prop.PT[, c("T", "P", "Cp")]
## ----Nalm, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na<sup>+</sup> (inapplicable: constant ω).", out.width=672, out.height=336, pngquant=pngquant----
var <- c("invTTheta2", "TXBorn")
Nalm <- EOSregress(Nadat, var, T.max = 600)
EOSplot(Nadat, coefficients = Nalm$coefficients, fun.legend = NULL)
EOSplot(Nadat, add = TRUE, lty = 3)
## ----Navars1------------------------------------------------------------------
var1 <- c("invTTheta2", "Cp_s_var")
omega.guess <- coef(Nalm)[3]
## ----Nawhile, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na<sup>+</sup> (variable ω).", out.width=672, out.height=336, pngquant=pngquant----
diff.omega <- 999
while(abs(diff.omega) > 1) {
Nalm1 <- EOSregress(Nadat, var1, omega.PrTr = tail(omega.guess, 1), Z = 1)
omega.guess <- c(omega.guess, coef(Nalm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(Nadat, coefficients = signif(coef(Nalm1), 6),
omega.PrTr = tail(omega.guess, 1), Z = 1)
convert(EOScoeffs("Na+", "Cp"), "J")
## ----NaVolume, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Volume of Na<sup>+</sup> (variable ω).", results="hide", message=FALSE, echo=FALSE, out.width=672, out.height=336, pngquant=pngquant----
T <- convert(seq(0, 600, 25), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
NaVdat <- prop.PT[, c("T", "P", "V")]
var1 <- c("invTTheta", "V_s_var")
omega.guess <- 1400000
diff.omega <- 999
while(abs(diff.omega) > 1) {
NaVlm1 <- EOSregress(NaVdat, var1,
omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1)
omega.guess <- c(omega.guess, coef(NaVlm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(NaVdat, coefficients = signif(coef(NaVlm1), 6),
omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1,
fun.legend = "bottomleft")
coefficients <- convert(EOScoeffs("Na+", "V", P = 1000), "J")
names(coefficients)[3] <- "V_s_var"
EOSplot(NaVdat, coefficients = coefficients, Z = 1, add = TRUE, lty = 2,
omega.PrTr = convert(coefficients["V_s_var"], "joules"))
## ----NaVolume, eval=FALSE-----------------------------------------------------
# T <- convert(seq(0, 600, 25), "K")
# P <- 1000
# prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
# NaVdat <- prop.PT[, c("T", "P", "V")]
# var1 <- c("invTTheta", "V_s_var")
# omega.guess <- 1400000
# diff.omega <- 999
# while(abs(diff.omega) > 1) {
# NaVlm1 <- EOSregress(NaVdat, var1,
# omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1)
# omega.guess <- c(omega.guess, coef(NaVlm1)[3])
# diff.omega <- tail(diff(omega.guess), 1)
# }
# EOSplot(NaVdat, coefficients = signif(coef(NaVlm1), 6),
# omega.PrTr = tail(convert(omega.guess, "joules"), 1), Z = 1,
# fun.legend = "bottomleft")
# coefficients <- convert(EOScoeffs("Na+", "V", P = 1000), "J")
# names(coefficients)[3] <- "V_s_var"
# EOSplot(NaVdat, coefficients = coefficients, Z = 1, add = TRUE, lty = 2,
# omega.PrTr = convert(coefficients["V_s_var"], "joules"))
## ----AS04, message=FALSE------------------------------------------------------
add.OBIGT("AS04")
## ----SiO2_2H2O, message=FALSE-------------------------------------------------
s_25C <- subcrt(c("SiO2", "H2O"), c(1, 2), T = 25)$out
s_near25 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(20, 30, length.out=50))$out
s_lowT <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 100, 10))$out
s_Psat <- subcrt(c("SiO2", "H2O"), c(1, 2))$out
s_P500 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 500)$out
s_P1000 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 1000)$out
## ----new_H4SiO4---------------------------------------------------------------
mod.OBIGT("calc-H4SiO4", formula = "H4SiO4", ref1 = "this_vignette",
date = as.character(Sys.Date()), G = s_25C$G, H = s_25C$H, S = s_25C$S,
Cp = s_25C$Cp, V = s_25C$V, z = 0)
## ----substuff-----------------------------------------------------------------
substuff <- rbind(s_near25, s_lowT, s_Psat, s_P500, s_P1000)
substuff$T <- convert(substuff$T, "K")
## ----Cp_H4SiO4, results="hide"------------------------------------------------
Cpdat <- substuff[, c("T", "P", "Cp")]
var <- c("invTTheta2", "TXBorn")
Cplm <- EOSregress(Cpdat, var)
Cpcoeffs <- Cplm$coefficients
mod.OBIGT("calc-H4SiO4", c1 = Cpcoeffs[1],
c2 = Cpcoeffs[2]/10000, omega = Cpcoeffs[3]/100000)
## ----V_H4SiO4_nonsolvation----------------------------------------------------
Vdat <- substuff[, c("T", "P", "V")]
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat$V <- V_non
## ----V_H4SiO4, results="hide"-------------------------------------------------
var <- c("invPPsi", "invTTheta", "invPPsiTTheta")
Vlm <- EOSregress(Vdat, var)
Vcoeffs <- convert(Vlm$coefficients, "joules")
mod.OBIGT("calc-H4SiO4", a1 = Vcoeffs[1]*10, a2 = Vcoeffs[2]/100,
a3 = Vcoeffs[3], a4 = Vcoeffs[4]/10000)
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----info_H4SiO4, message=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------
info(info(c("calc-H4SiO4", "H4SiO4")))
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----subcrt_H4SiO4, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=FALSE, results="hide", message=FALSE, out.width="100%", cache=TRUE, fig.cap="Comparison of H<sub>4</sub>SiO<sub>4</sub> pseudospecies.", pngquant=pngquant----
s1 <- subcrt(c("calc-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
plot(s1$out$T, s1$out$G, type = "l", ylim = c(-500, 2000),
xlab = axis.label("T"), ylab = axis.label("DG0"))
s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
lines(s2$out$T, s2$out$G, lty = 2)
abline(h = 0, lty = 3)
legend("topright", legend = c("calc-H4SiO4 (this vignette)",
"H4SiO4 (Stef\u00e1nsson, 2001)"), lty = c(1, 2), bty = "n")
text(225, 1000, describe.reaction(s1$reaction))
## ----subcrt_H4SiO4, eval=FALSE------------------------------------------------
# s1 <- subcrt(c("calc-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
# plot(s1$out$T, s1$out$G, type = "l", ylim = c(-500, 2000),
# xlab = axis.label("T"), ylab = axis.label("DG0"))
# s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
# lines(s2$out$T, s2$out$G, lty = 2)
# abline(h = 0, lty = 3)
# legend("topright", legend = c("calc-H4SiO4 (this vignette)",
# "H4SiO4 (Stef\u00e1nsson, 2001)"), lty = c(1, 2), bty = "n")
# text(225, 1000, describe.reaction(s1$reaction))
## ----activity_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=TRUE, results="hide", message=FALSE, out.width="100%", cache=TRUE, fig.cap="Activity diagram for K<sub>2</sub>O-Al<sub>2</sub>O<sub>3</sub>-SiO<sub>2</sub>-H<sub>2</sub>O.", pngquant=pngquant----
basis(c("Al+3", "H4SiO4", "K+", "H2O", "H+", "O2"))
species(c("gibbsite", "muscovite", "kaolinite", "pyrophyllite", "K-feldspar"))
a <- affinity(H4SiO4 = c(-8, 0, 300), `K+` = c(-1, 8, 300))
diagram(a, ylab = ratlab("K+"), fill = "terrain", yline = 1.7)
legend("bottomleft", describe.property(c("T", "P"), c(25, 1)), bty = "n")
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.