# demos/buffer.R
# calculate buffered activities of basis species using two methods
# reproduces Figure 6 of Schulte and Shock, 1995 (doi:10.1007/BF01581580)

b.species <- c("Fe", "CO2", "H2O", "N2", "H2", "H2S", "SiO2")
b.state <- c("cr", "gas", "liq", "gas", "gas", "aq", "aq")
b.logact <- c(0, 1, 0, 0, 0, 0, 0)
basis(b.species, b.state, b.logact)
xlim <- c(0, 350)
thermo.plot.new(xlim=xlim, ylim=c(-4, 4), xlab=axis.label("T"), ylab=axis.label("H2"))
# method 1: in affinity(), assign name of buffer to basis species
bufferline <- function(buffer, ixlab) {
  basis("H2", buffer)
  a <- affinity(T=xlim, P=300, return.buffer=TRUE, exceed.Ttr=TRUE)
  lines(a$vals[[1]], a$H2, col=3, lwd=2)
  text(a$vals[[1]][ixlab], a$H2[ixlab] + 0.2, buffer)
}
bufferline("FeFeO", 20)
bufferline("QFM", 38)
bufferline("PPM", 102)
bufferline("HM", 51)
# method 2: in diagram(), use the `what` argument
basis("H2", 0)
for(logact in c(-6, -10, -15)) {
  species(c("formaldehyde", "HCN"), logact)
  a <- affinity(T=xlim, P=300)
  d <- diagram(a, what="H2", lty=c(3, 2), add=TRUE)
  text(a$vals[[1]][13], mean(sapply(d$plotvals, c)[13, ]), logact)
}
# add legends and title
legend("topright", legend = c("minerals", "formaldehyde", "HCN"),
  lty=c(1, 3, 2), col=c(3, 1, 1), bg="white", cex=0.9)
legend("bottomright", legend = c(describe.property("P", 300),
  describe.basis(ibasis=c(2,4))), bg="white", cex=0.9)
title(main=paste("Mineral buffers and activities of aqueous species",
                 "(Schulte and Shock, 1995)", sep="\n"), cex.main=0.9)
