[CHNOSZ-commits] r770 - in pkg/CHNOSZ: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 2 11:06:04 CET 2023
Author: jedick
Date: 2023-03-02 11:06:03 +0100 (Thu, 02 Mar 2023)
New Revision: 770
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/DEW.R
pkg/CHNOSZ/R/IAPWS95.R
pkg/CHNOSZ/R/add.OBIGT.R
pkg/CHNOSZ/R/add.protein.R
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/ionize.aa.R
pkg/CHNOSZ/R/makeup.R
pkg/CHNOSZ/R/mix.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/palply.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/retrieve.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/swap.basis.R
pkg/CHNOSZ/R/taxonomy.R
pkg/CHNOSZ/R/thermo.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/util.array.R
pkg/CHNOSZ/R/util.character.R
pkg/CHNOSZ/R/util.fasta.R
pkg/CHNOSZ/R/util.formula.R
pkg/CHNOSZ/R/util.legend.R
pkg/CHNOSZ/R/util.list.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/R/util.protein.R
pkg/CHNOSZ/R/util.seq.R
pkg/CHNOSZ/R/util.units.R
pkg/CHNOSZ/R/util.water.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/R/zzz.R
Log:
Capitalize comments
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/DESCRIPTION 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,6 +1,6 @@
Date: 2023-03-02
Package: CHNOSZ
-Version: 1.9.9-61
+Version: 1.9.9-62
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/DEW.R
===================================================================
--- pkg/CHNOSZ/R/DEW.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/DEW.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -15,7 +15,7 @@
# 'pressure - The pressure to calculate the density of water at, in bars
# 'temperature - The temperature to calculate the density of water at, in Celsius
# 'error - The density returned will calculate a pressure which differs from the input pressure by the value of "error" or less.
-# the default value of error is taken from equations in DEW spreadsheet (table "Calculations")
+# The default value of error is taken from equations in DEW spreadsheet (table "Calculations")
calculateDensity <- function(pressure, temperature, error = 0.01) {
myfunction <- function(pressure, temperature) {
minGuess <- 1E-05
@@ -46,10 +46,10 @@
}
calculateDensity
}
- # make input pressure and temperature the same length
+ # Make input pressure and temperature the same length
if(length(pressure) < length(temperature)) pressure <- rep(pressure, length.out=length(temperature))
if(length(temperature) < length(pressure)) temperature <- rep(temperature, length.out=length(pressure))
- # use a loop to process vectorized input
+ # Use a loop to process vectorized input
sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
}
@@ -87,10 +87,10 @@
}
GAtOneKb + integral
}
- # make input pressure and temperature the same length
+ # Make input pressure and temperature the same length
if(length(pressure) < length(temperature)) pressure <- rep(pressure, length.out=length(temperature))
if(length(temperature) < length(pressure)) temperature <- rep(temperature, length.out=length(pressure))
- # use a loop to process vectorized input
+ # Use a loop to process vectorized input
sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
}
@@ -125,7 +125,7 @@
depsdrho * drhodP / eps ^2
}
-### unexported functions ###
+### Unexported functions ###
# 'Returns the pressure of water corresponding to the input density and temperature, in units of bars.
# 'density - The density to use in finding a pressure, in g/cm^3
@@ -200,7 +200,7 @@
A * exp(B) * density ^ (A - 1)
}
-### testing functions ###
+### Testing functions ###
# These unexported functions are included for testing purposes only.
# In CHNOSZ, the g function and omega(P,T) are calculated via hkf().
@@ -241,7 +241,7 @@
(-1.504956E-10 * (1000 - P)^3 + 5.017997E-14 * (1000 - P)^4)
f[P > 1000 | T < 155 | T > 355] <- 0
g <- a_g * (1 - density)^b_g - f
- # use g = 0 for density >= 1
+ # Use g = 0 for density >= 1
g[density >= 1] <- 0
g
}
Modified: pkg/CHNOSZ/R/IAPWS95.R
===================================================================
--- pkg/CHNOSZ/R/IAPWS95.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/IAPWS95.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,32 +1,29 @@
-# functions for properties of water using
-# the IAPWS-95 formulation (Wagner and Pruss, 2002)
+# Calculate properties of water using the IAPWS-95 formulation (Wagner and Pruss, 2002)
IAPWS95 <- function(property,T=298.15,rho=1000) {
- ## the IAPWS-95 formulation for ordinary water substance
- ## Wagner and Pruss, 2002
property <- tolower(property)
- # triple point
+ # Triple point
T.triple <- 273.16 # K
P.triple <- 611.657 # Pa
rho.triple.liquid <- 999.793
rho.triple.vapor <- 0.00485458
- # normal boiling point
+ # Normal boiling point
T.boiling <- 373.124
P.boiling <- 0.101325
rho.boiling.liquid <- 958.367
rho.boiling.vapor <- 0.597657
- # critical point constants
+ # Critical point constants
T.critical <- 647.096 # K
rho.critical <- 322 # kg m-3
- # specific and molar gas constants
+ # Specific and molar gas constants
R <- 0.46151805 # kJ kg-1 K-1
# R.M <- 8.314472 # J mol-1 K-1
- # molar mass
+ # Molar mass
M <- 18.015268 # g mol-1
- ## define functions idealgas and residual, supplying arguments delta and tau
+ ## Define functions idealgas and residual, supplying arguments delta and tau
idealgas <- function(p) IAPWS95.idealgas(p, delta, tau)
residual <- function(p) IAPWS95.residual(p, delta, tau)
- ## relation of thermodynamic properties to Helmholtz free energy
+ ## Relation of thermodynamic properties to Helmholtz free energy
a <- function() {
x <- idealgas('phi')+residual('phi')
return(x*R*T)
@@ -62,7 +59,7 @@
(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta'))
return(x*R)
}
-# 20090420 speed of sound calculation is incomplete
+# 20090420 Speed of sound calculation is incomplete
# (delta.liquid and drhos.dT not visible)
# cs <- function() {
# x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
@@ -83,7 +80,7 @@
(idealgas('phi.tau.tau')+residual('phi.tau.tau'))*(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) )
return(x/(R*rho))
}
- ## run the calculations
+ ## Run the calculations
ww <- NULL
my.T <- T
my.rho <- rho
@@ -104,7 +101,7 @@
return(ww)
}
-### unexported functions ###
+### Unexported functions ###
# IAPWS95.idealgas and IAPWS95.residual are supporting functions to IAPWS95 for calculating
# the ideal-gas and residual parts in the IAPWS-95 formulation.
@@ -112,8 +109,8 @@
# to calculate the specific dimensionless Helmholtz free energy (phi) or one of its derivatives.
IAPWS95.idealgas <- function(p, delta, tau) {
- ## the ideal gas part in the IAPWS-95 formulation
- # from Table 6.1 of Wagner and Pruss, 2002
+ ## The ideal gas part in the IAPWS-95 formulation
+ # From Table 6.1 of Wagner and Pruss, 2002
n <- c( -8.32044648201, 6.6832105268, 3.00632, 0.012436,
0.97315, 1.27950, 0.96956, 0.24873 )
gamma <- c( NA, NA, NA, 1.28728967,
@@ -121,7 +118,7 @@
# Equation 6.5
phi <- function() log(delta) + n[1] + n[2]*tau + n[3]*log(tau) +
sum( n[4:8] * log(1-exp(-gamma[4:8]*tau)) )
- # derivatives from Table 6.4
+ # Derivatives from Table 6.4
phi.delta <- function() 1/delta+0+0+0+0
phi.delta.delta <- function() -1/delta^2+0+0+0+0
phi.tau <- function() 0+0+n[2]+n[3]/tau+sum(n[4:8]*gamma[4:8]*((1-exp(-gamma[4:8]*tau))^-1-1))
@@ -132,8 +129,8 @@
}
IAPWS95.residual <- function(p, delta, tau) {
- ## the residual part in the IAPWS-95 formulation
- # from Table 6.2 of Wagner and Pruss, 2002
+ ## The residual part in the IAPWS-95 formulation
+ # From Table 6.2 of Wagner and Pruss, 2002
c <- c(rep(NA,7),rep(1,15),rep(2,20),rep(3,4),4,rep(6,4),rep(NA,5))
d <- c(1,1,1,2,2,3,4,1,1,1,2,2,3,4,
4,5,7,9,10,11,13,15,1,2,2,2,3,4,
@@ -177,7 +174,7 @@
i2 <- 8:51
i3 <- 52:54
i4 <- 55:56
- # deriviatives of distance function
+ # Deriviatives of distance function
Delta <- function(i) { Theta(i)^2 + B[i] * ((delta-1)^2)^a[i] }
Theta <- function(i) { (1-tau) + A[i] * ((delta-1)^2)^(1/(2*beta[i])) }
Psi <- function(i) { exp ( -C[i]*(delta-1)^2 - D[i]*(tau-1)^2 ) }
@@ -194,13 +191,13 @@
4*B[i]*a[i]*(a[i]-1)*((delta-1)^2)^(a[i]-2) + 2*A[i]^2*(1/beta[i])^2 *
(((delta-1)^2)^(1/(2*B[i])-1))^2 + A[i]*Theta(i)*4/beta[i]*(1/(2*B[i])-1) *
((delta-1)^2)^(1/(2*beta[i])-2) ) }
- # derivatives of exponential function
+ # Derivatives of exponential function
dPsi.ddelta <- function(i) { -2*C[i]*(delta-1)*Psi(i) }
d2Psi.ddelta2 <- function(i) { ( 2*C[i]*(delta-1)^2 - 1 ) * 2*C[i]*Psi(i) }
dPsi.dtau <- function(i) { -2*D[i]*(tau-1)*Psi(i) }
d2Psi.dtau2 <- function(i) { (2*D[i]*(tau-1)^2 - 1) * 2*D[i]*Psi(i) }
d2Psi.ddelta.dtau <- function(i) { 4*C[i]*D[i]*(delta-1)*(tau-1)*Psi(i) }
- # dimensionless Helmholtz free energy and derivatives
+ # Dimensionless Helmholtz free energy and derivatives
phi <- function() {
sum(n[i1]*delta^d[i1]*tau^t[i1]) +
sum(n[i2]*delta^d[i2]*tau^t[i2]*exp(-delta^c[i2])) +
Modified: pkg/CHNOSZ/R/add.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/add.OBIGT.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/add.OBIGT.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,5 +1,5 @@
# CHNOSZ/add.OBIGT.R
-# add or change entries in the thermodynamic database
+# Add or change entries in the thermodynamic database
## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("info.R")
Modified: pkg/CHNOSZ/R/add.protein.R
===================================================================
--- pkg/CHNOSZ/R/add.protein.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/add.protein.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -4,7 +4,7 @@
# Count numbers of amino acids in a sequence
seq2aa <- function(sequence, protein = NA) {
- # $emove newlines and whitespace
+ # Remove newlines and whitespace
sequence <- gsub("\\s", "", gsub("[\r\n]", "", sequence))
# Make a data frame from counting the amino acids in the sequence
caa <- count.aa(sequence)
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/affinity.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,7 +1,7 @@
# CHNOSZ/affinity.R
-# calculate affinities of formation reactions
+# Calculate affinities of formation reactions
-## if this file is interactively sourced, the following are also needed to provide unexported functions:
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("util.affinity.R")
#source("util.units.R")
#source("util.character.R")
@@ -24,19 +24,19 @@
# sout: provide a previously calculated output from subcrt
# iprotein: build these proteins from residues (speed optimization)
- # history: 20061027 jmd version 1
+ # History: 20061027 jmd version 1
# this is where energy.args() used to sit
# this is where energy() used to sit
- # argument recall 20190117
- # if the first argument is the result from a previous affinity() calculation,
+ # Argument recall 20190117
+ # If the first argument is the result from a previous affinity() calculation,
# just update the remaining arguments
args.orig <- list(...)
- # we can only do anything with at least one argument
+ # We can only do anything with at least one argument
if(length(args.orig) > 0) {
if(identical(args.orig[[1]][1], list(fun="affinity"))) {
aargs <- args.orig[[1]]$args
- # we can only update arguments given after the first argument
+ # We can only update arguments given after the first argument
if(length(args.orig) > 1) {
for(i in 2:length(args.orig)) {
if(names(args.orig)[i] %in% names(aargs)) aargs[[names(args.orig)[i]]] <- args.orig[[i]]
@@ -47,17 +47,17 @@
}
}
- # the argument list
+ # The argument list
args <- energy.args(args.orig, transect = transect)
args <- c(args, list(sout=sout, exceed.Ttr=exceed.Ttr, exceed.rhomin=exceed.rhomin))
- # the species we're given
+ # The species we're given
thermo <- get("thermo", CHNOSZ)
mybasis <- thermo$basis
myspecies <- thermo$species
if(!is.null(property)) {
- # the user just wants an energy property
+ # The user just wants an energy property
buffer <- FALSE
args$what <- property
out <- do.call("energy",args)
@@ -65,21 +65,18 @@
sout <- out$sout
} else {
- # affinity calculations
+ # Affinity calculations
property <- args$what
# iprotein stuff
- # note that this affinities of the residues are subject
- # to ionization calculations in energy() so no explicit accounting
- # is needed here
+ # Note that affinities of the residues are modified by ionization calculations in energy(), not here
if(!is.null(iprotein)) {
- # check all proteins are available
+ # Check all proteins are available
if(any(is.na(iprotein))) stop("`iprotein` has some NA values")
if(!all(iprotein %in% 1:nrow(thermo$protein))) stop("some value(s) of `iprotein` are not rownumbers of thermo()$protein")
- # add protein residues to the species list
+ # Add protein residues to the species list
resnames <- c("H2O",aminoacids(3))
- # residue activities set to zero;
- # account for protein activities later
+ # Residue activities set to zero; account for protein activities later
resprot <- paste(resnames,"RESIDUE",sep="_")
species(resprot, 0)
thermo <- get("thermo", CHNOSZ)
@@ -86,10 +83,10 @@
ires <- match(resprot, thermo$species$name)
}
- # buffer stuff
+ # Buffer stuff
buffer <- FALSE
hasbuffer <- !can.be.numeric(mybasis$logact)
- # however, variables specified in the arguments take precedence over the buffer 20171013
+ # However, variables specified in the arguments take precedence over the buffer 20171013
isnotvar <- !rownames(mybasis) %in% args$vars
ibufbasis <- which(hasbuffer & isnotvar)
if(!is.null(mybasis) & length(ibufbasis) > 0) {
@@ -102,11 +99,11 @@
for(i in 1:length(is.buffer)) is.buff <- c(is.buff,as.numeric(is.buffer[[i]]))
is.only.buffer <- is.buff[!is.buff %in% is.species]
buffers <- names(is.buffer)
- # reorder the buffers according to thermo$buffer
+ # Reorder the buffers according to thermo$buffer
buffers <- buffers[order(match(buffers,thermo$buffer$name))]
}
- # here we call 'energy'
+ # Here we call 'energy'
aa <- do.call("energy",args)
a <- aa$a
sout <- aa$sout
@@ -113,7 +110,7 @@
if(return.sout) return(sout)
- # more buffer stuff
+ # More buffer stuff
if(buffer) {
args$what <- "logact.basis"
args$sout <- sout
@@ -121,7 +118,7 @@
ibasis.new <- numeric()
for(k in 1:length(buffers)) {
ibasis <- which(as.character(mybasis$logact)==buffers[k])
- # calculate the logKs from the affinities
+ # Calculate the logKs from the affinities
logK <- a
for(i in 1:length(logK)) {
logK[[i]] <- logK[[i]] + thermo$species$logact[i]
@@ -132,12 +129,12 @@
lbn <- buffer(logK=logK,ibasis=ibasis,logact.basis=logact.basis.new,
is.buffer=as.numeric(is.buffer[[which(names(is.buffer)==buffers[k])]]),balance=balance)
for(j in 1:length(logact.basis.new)) if(j %in% ibasis) logact.basis.new[[j]] <- lbn[[2]][[j]]
- # calculation of the buffered activities' effect on chemical affinities
+ # Calculation of the buffered activities' effect on chemical affinities
is.only.buffer.new <- is.only.buffer[is.only.buffer %in% is.buffer[[k]]]
for(i in 1:length(a)) {
if(i %in% is.only.buffer.new) next
for(j in 1:nrow(mybasis)) {
- # let's only do this for the basis species specified by the user
+ # Let's only do this for the basis species specified by the user
# even if others could be buffered
if(!j %in% ibufbasis) next
if(!j %in% ibasis) next
@@ -153,7 +150,7 @@
}
species(is.only.buffer,delete=TRUE)
if(length(is.only.buffer) > 0) a <- a[-is.only.buffer]
- # to return the activities of buffered basis species
+ # To return the activities of buffered basis species
tb <- logact.basis.new[unique(ibasis.new)]
if(!is.null(ncol(tb[[1]]))) {
nd <- sum(dim(tb[[1]]) > 1)
@@ -170,38 +167,36 @@
}
}
- # more iprotein stuff
+ # More iprotein stuff
if(!is.null(iprotein)) {
- # 20090331 fast protein calculations
- # function to calculate affinity of formation reactions
- # from those of residues
+ # Fast protein calculations 20090331
+ # Function to calculate affinity of formation reactions from those of residues
loga.protein <- rep(loga.protein,length.out=length(iprotein))
protein.fun <- function(ip) {
tpext <- as.numeric(thermo$protein[iprotein[ip],5:25])
return(Reduce("+", pprod(a[ires],tpext)) - loga.protein[ip])
}
- # use another level of indexing to let the function
- # report on its progress
+ # Use another level of indexing to let the function report on its progress
jprotein <- 1:length(iprotein)
protein.affinity <- palply("", jprotein, protein.fun)
- ## update the species list
- # we use negative values for ispecies to denote that
+ ## Update the species list
+ # We use negative values for ispecies to denote that
# they index thermo$protein and not thermo$species
ispecies <- -iprotein
- # the current species list, containing the residues
+ # The current species list, containing the residues
resspecies <- thermo$species
- # now we can delete the residues from the species list
+ # Now we can delete the residues from the species list
species(ires,delete=TRUE)
- # state and protein names
+ # State and protein names
state <- resspecies$state[1]
name <- paste(thermo$protein$protein[iprotein],thermo$protein$organism[iprotein],sep="_")
- # the numbers of basis species in formation reactions of the proteins
+ # The numbers of basis species in formation reactions of the proteins
protbasis <- t(t((resspecies[ires,1:nrow(mybasis)])) %*% t((thermo$protein[iprotein,5:25])))
- # put them together
+ # Put them together
protspecies <- cbind(protbasis,data.frame(ispecies=ispecies,logact=loga.protein,state=state,name=name))
myspecies <- rbind(myspecies,protspecies)
rownames(myspecies) <- 1:nrow(myspecies)
- ## update the affinity values
+ ## Update the affinity values
names(protein.affinity) <- ispecies
a <- c(a,protein.affinity)
a <- a[-ires]
@@ -209,14 +204,14 @@
}
- # put together return values
- # constant T and P, in Kelvin and bar
+ # Put together return values
+ # Constant T and P, in Kelvin and bar
T <- args$T
P <- args$P
- # the variable names and values
+ # The variable names and values
vars <- args$vars
vals <- args$vals
- # if T or P is variable it's not constant;
+ # If T or P is variable it's not constant;
# and set it to user's units
iT <- match("T", vars)
if(!is.na(iT)) {
@@ -228,11 +223,11 @@
P <- numeric()
vals[[iP]] <- outvert(vals[[iP]], "bar")
}
- # get Eh
+ # Get Eh
iEh <- match("Eh", names(args.orig))
if(!is.na(iEh)) {
vars[[iEh]] <- "Eh"
- # we have to reconstruct the values used because
+ # We have to reconstruct the values used because
# they got converted to log_a(e-) at an unknown temperature
Eharg <- args.orig[[iEh]]
if(length(Eharg) > 3) Ehvals <- Eharg
@@ -240,7 +235,7 @@
else if(length(Eharg) == 2) Ehvals <- seq(Eharg[1], Eharg[2], length.out=256)
vals[[iEh]] <- Ehvals
}
- # get pe and pH
+ # Get pe and pH
ipe <- match("pe", names(args.orig))
if(!is.na(ipe)) {
ie <- match("e-", names(args$lims))
@@ -253,15 +248,15 @@
vars[[iH]] <- "pH"
vals[[iH]] <- -args$vals[[iH]]
}
- # use the variable names for the vals list 20190415
+ # Use the variable names for the vals list 20190415
names(vals) <- vars
- # content of return value depends on buffer request
+ # Content of return value depends on buffer request
if(return.buffer) return(c(tb, list(vars=vars, vals=vals)))
- # for argument recall, include all arguments (except sout) in output 20190117
+ # For argument recall, include all arguments (except sout) in output 20190117
allargs <- c(args.orig, list(property=property, exceed.Ttr=exceed.Ttr, exceed.rhomin=exceed.rhomin,
return.buffer=return.buffer, balance=balance, iprotein=iprotein, loga.protein=loga.protein))
- # add IS value only if it given as an argument 20171101
+ # Add IS value only if it given as an argument 20171101
# (even if its value is 0, the presence of IS will trigger diagram() to use "m" instead of "a" in axis labels)
iIS <- match("IS", names(args.orig))
if(!is.na(iIS)) a <- list(fun="affinity", args=allargs, sout=sout, property=property,
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/buffer.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -56,14 +56,14 @@
return(invisible(thermo$buffer[thermo$buffer$name %in% name,]))
}
-### unexported functions ###
+### Unexported functions ###
buffer <- function(logK=NULL,ibasis=NULL,logact.basis=NULL,is.buffer=NULL,balance='PBB') {
thermo <- get("thermo", CHNOSZ)
- # if logK is NULL load the buffer species
- # otherwise perform buffer calculations.
+ # If logK is NULL load the buffer species,
+ # otherwise perform buffer calculations
if(is.null(logK)) {
- # load the buffer species
+ # Load the buffer species
buffers <- unique(as.character(thermo$basis$logact)[!can.be.numeric(as.character(thermo$basis$logact))])
ispecies.new <- list()
for(k in 1:length(buffers)) {
@@ -77,7 +77,7 @@
ispecies <- c(ispecies, species(species, state, index.return=TRUE, add = TRUE))
}
ispecies.new <- c(ispecies.new,list(ispecies))
- # make sure to set the activities
+ # Make sure to set the activities
species(ispecies,thermo$buffer$logact[ib])
}
names(ispecies.new) <- buffers
@@ -84,9 +84,9 @@
return(ispecies.new)
}
- # sometimes (e.g. PPM) the buffer species are identified multiple
- # times, causing problems for square matrices and such.
- # make them appear singly.
+ # Sometimes (e.g. PPM) the buffer species are identified multiple
+ # times, causing problems for square matrices
+ # Make the species appear singly
is.buffer <- unique(is.buffer)
bufbasis <- species.basis(thermo$species$ispecies[is.buffer])
bufname <- thermo$basis$logact[ibasis[1]]
@@ -99,29 +99,29 @@
bufbasis <- cbind(data.frame(product=nb),bufbasis)
} else {
basisnames <- c('PBB',basisnames)
- # prepend a PBB column to bufbasis and inc. ibasis by 1
+ # Prepend a PBB column to bufbasis and increment ibasis by 1
nb <- as.numeric(protein.length(thermo$species$name[is.buffer]))
bufbasis <- cbind(data.frame(PBB=nb),bufbasis)
}
ibasis <- ibasis + 1
- # make logact.basis long enough
+ # Make logact.basis long enough
tl <- length(logact.basis)
logact.basis[[tl+1]] <- logact.basis[[tl]]
- # rotate the entries so that the new one is first
+ # Rotate the entries so that the new one is first
ilb <- c(tl+1,1:tl)
logact.basis <- logact.basis[ilb]
}
- # say hello
- #cat(paste("buffer: '",bufname,"', of ",length(is.buffer),' species, ',length(ibasis),' activity(s) requested.\n',sep=''))
+ # Say hello
+ #message(paste("buffer: '",bufname,"', of ",length(is.buffer),' species, ',length(ibasis),' activity(s) requested.',sep=''))
ibasisrequested <- ibasis
- # check and maybe add to the number of buffered activities
+ # Check and maybe add to the number of buffered activities
ibasisadded <- numeric()
if( (length(ibasis)+1) != length(is.buffer) & length(is.buffer) > 1) {
- # try to add buffered activities the user didn't specify
+ # Try to add buffered activities the user didn't specify
# (e.g. H2S in PPM buffer if only H2 was requested)
for(i in 1:(length(is.buffer)-(length(ibasis)+1))) {
newbasis <- NULL
- # we want to avoid any basis species that might be used as the conservant
+ # We want to avoid any basis species that might be used as the conservant
# look for additional activities to buffer ... do columns in reverse
for(j in ncol(bufbasis):1) {
if(j %in% ibasis) next
@@ -138,17 +138,17 @@
}
}
}
- # and the leftovers
+ # And the leftovers
#xx <- as.data.frame(bufbasis[,-ibasis])
- # the final buffered activity: the would-be conserved component
+ # The final buffered activity: the would-be conserved component
newbasis <- NULL
if(length(is.buffer) > 1) {
- # first try to get one that is present in all species
+ # First try to get one that is present in all species
for(i in ncol(bufbasis):1) {
if(i %in% ibasis) next
if(!TRUE %in% (bufbasis[,i]==0)) newbasis <- i
}
- # or look for one that is present at all
+ # Or look for one that is present at all
if(is.null(newbasis)) for(i in ncol(bufbasis):1) {
if(i %in% ibasis) next
if(FALSE %in% (bufbasis[,i]==0)) newbasis <- i
@@ -155,7 +155,7 @@
}
if(!is.null(newbasis)) {
ibasis <- c(ibasis,newbasis)
- #cat(paste('buffer: the conserved activity is ',basisnames[newbasis],'.\n',sep=''))
+ #message(paste('buffer: the conserved activity is ',basisnames[newbasis],'.',sep=''))
#thermo$basis$logact[newbasis] <<- thermo$basis$logact[ibasis[1]]
}
else stop('no conserved activity found in your buffer (not enough basis species?)!')
@@ -164,15 +164,14 @@
reqtext <- paste(c2s(basisnames[ibasisrequested]),' (active)',sep='')
if(length(ibasisadded)==0) addtext <- '' else addtext <- paste(', ',c2s(basisnames[ibasisadded]),sep='')
message(paste("buffer: '", bufname, "' for activity of ", reqtext, addtext, context, sep = ""))
- #print(bufbasis)
- # there could still be stuff here (over-defined system?)
+ # There could still be stuff here (over-defined system?)
xx <- bufbasis[,-ibasis,drop=FALSE]
- # for the case when all activities are buffered
+ # For the case when all activities are buffered
if(ncol(xx)==0) xx <- data.frame(xx=0)
- # our stoichiometric matrix - should be square
+ # Our stoichiometric matrix - should be square
A <- as.data.frame(bufbasis[,ibasis])
- # determine conservation coefficients
- # values for the known vector
+ # Determine conservation coefficients
+ # Values for the known vector
B <- list()
for(i in 1:length(is.buffer)) {
b <- -logK[[is.buffer[i]]] + thermo$species$logact[is.buffer[i]]
@@ -179,7 +178,7 @@
if(ncol(xx) > 0) {
if(is.list(xx)) xxx <- xx[[1]] else xxx <- xx
if(ncol(xx)==1 & identical(as.numeric(xxx),0)) {
- # do nothing
+ # Do nothing
} else {
for(j in 1:ncol(xx)) {
bs <- xx[i,j] * logact.basis[[match(colnames(xx)[j],basisnames)]]
@@ -192,7 +191,7 @@
# Force this to be matrix so indexing works at a single point (no variables defined in affinity()) 20201102
B[[i]] <- as.matrix(b)
}
- # a place to put the results
+ # A place to put the results
X <- rep(B[1],length(ibasis))
for(i in 1:nrow(B[[1]])) {
for(j in 1:ncol(B[[1]])) {
@@ -205,7 +204,7 @@
X[[k]][i,j] <- t[k]
}
}
- # store results
+ # Store results
for(i in 1:length(ibasis)) {
if(ncol(X[[i]])==1) X[[i]] <- as.numeric(X[[i]])
else if(nrow(X[[i]])==1) X[[i]] <- as.matrix(X[[i]],nrow=1)
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/cgl.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -95,10 +95,10 @@
### Unexported function ###
-# calculate GHS and V corrections for quartz and coesite 20170929
+# Calculate GHS and V corrections for quartz and coesite 20170929
# (these are the only mineral phases for which SUPCRT92 uses an inconstant volume)
quartz_coesite <- function(PAR, T, P) {
- # the corrections are 0 for anything other than quartz and coesite
+ # The corrections are 0 for anything other than quartz and coesite
if(!PAR$name %in% c("quartz", "coesite")) return(list(G=0, H=0, S=0, V=0))
ncond <- max(c(length(T), length(P)))
# Tr, Pr and TtPr (transition temperature at Pr)
@@ -105,7 +105,7 @@
Pr <- 1 # bar
Tr <- 298.15 # K
TtPr <- 848 # K
- # constants from SUP92D.f
+ # Constants from SUP92D.f
aa <- 549.824
ba <- 0.65995
ca <- -0.4973e-4
@@ -112,12 +112,12 @@
VPtTta <- 23.348
VPrTtb <- 23.72
Stran <- convert(0.342, "J")
- # constants from REAC92D.f
+ # Constants from REAC92D.f
VPrTra <- 22.688 # VPrTr(a-quartz)
Vdiff <- 2.047 # VPrTr(a-quartz) - VPrTr(coesite)
k <- 38.5 # dPdTtr(a/b-quartz)
#k <- 38.45834 # calculated in CHNOSZ: dPdTtr(info("quartz"))
- # code adapted from REAC92D.f
+ # Code adapted from REAC92D.f
qphase <- gsub("cr", "", PAR$state)
if(qphase == 2) {
Pstar <- P
@@ -143,7 +143,6 @@
k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/k)))
SVterm <- cm3bar_to_J * (-k * (ba + aa * ca * k) *
log((aa + P/k) / (aa + Pstar/k)) + ca * k * (P - Pstar)) - Sstar
- # note the minus sign on "SVterm" in order that intdVdTdP has the correct sign
+ # Note the minus sign on "SVterm" in order that intdVdTdP has the correct sign
list(intVdP=GVterm, intdVdTdP=-SVterm, V=V)
}
-
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/diagram.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,11 +1,11 @@
# CHNOSZ/diagram.R
-# plot equilibrium chemical activity and predominance diagrams
+# Plot equilibrium chemical activity and predominance diagrams
# 20061023 jmd v1
# 20120927 work with output from either equil() or affinity(),
# gather plotvals independently of plot parameters (including nd),
# single return statement
-## if this file is interactively sourced, the following are also needed to provide unexported functions:
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("equilibrate.R")
#source("util.plot.R")
#source("util.character.R")
@@ -12,34 +12,34 @@
#source("util.misc.R")
diagram <- function(
- # species affinities or activities
+ # Species affinities or activities
eout,
- # type of plot
+ # Type of plot
type="auto", alpha=FALSE, normalize=FALSE, as.residue=FALSE, balance=NULL,
groups=as.list(1:length(eout$values)), xrange=NULL,
- # figure size and sides for axis tick marks
+ # Figure size and sides for axis tick marks
mar=NULL, yline=par("mgp")[1]+0.3, side=1:4,
- # axis limits and labels
+ # Axis limits and labels
ylog=TRUE, xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL,
- # sizes
+ # Sizes
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
- # line styles
+ # Line styles
lty=NULL, lty.cr=NULL, lty.aq=NULL, lwd=par("lwd"), dotted=NULL,
spline.method=NULL, contour.method="edge", levels=NULL,
- # colors
+ # Colors
col=par("col"), col.names=par("col"), fill=NULL,
fill.NA="gray80", limit.water=NULL,
- # field and line labels
+ # Field and line labels
names=NULL, format.names=TRUE, bold=FALSE, italic=FALSE,
font=par("font"), family=par("family"), adj=0.5, dx=0, dy=0, srt=0,
min.area=0,
- # title and legend
+ # Title and legend
main=NULL, legend.x=NA,
- # plotting controls
+ # Plotting controls
add=FALSE, plot.it=TRUE, tplot=TRUE, ...
) {
- ### argument handling ###
+ ### Argument handling ###
## Check that eout is a valid object
efun <- eout$fun
@@ -68,7 +68,7 @@
if(type %in% c("auto", "saturation")) {
if(!"loga.equil" %in% names(eout)) {
eout.is.aout <- TRUE
- # get the balancing coefficients
+ # Get the balancing coefficients
if(type=="auto") {
bal <- balance(eout, balance)
n.balance <- bal$n.balance
@@ -237,7 +237,7 @@
if(!wl$swapped) {
# The actual calculation
iNA <- eout$vals[[2]] < ymin | eout$vals[[2]] > ymax
- # assign NA to the predominance matrix
+ # Assign NA to the predominance matrix
H2O.predominant[i, iNA] <- NA
} else {
# As above, but x- and y-axes are swapped
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2023-03-02 09:02:20 UTC (rev 769)
+++ pkg/CHNOSZ/R/equilibrate.R 2023-03-02 10:06:03 UTC (rev 770)
@@ -1,8 +1,8 @@
# CHNOSZ/equilibrate.R
-# functions to calculation logarithm of activity
+# Functions to calculation logarithm of activity
# of species in (metastable) equilibrium
-## if this file is interactively sourced, the following are also needed to provide unexported functions:
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("util.misc.R")
#source("util.character.R")
@@ -9,48 +9,48 @@
equilibrate <- function(aout, balance=NULL, loga.balance=NULL,
ispecies = !grepl("cr", aout$species$state), normalize=FALSE, as.residue=FALSE,
method=c("boltzmann", "reaction"), tol=.Machine$double.eps^0.25) {
- ### calculate equilibrium activities of species from the affinities
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 770
More information about the CHNOSZ-commits
mailing list