Skip to content

Commit

Permalink
Add spaces around equals signs
Browse files Browse the repository at this point in the history
git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@793 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
  • Loading branch information
jedick committed Jun 24, 2023
1 parent 1978551 commit aa13309
Show file tree
Hide file tree
Showing 48 changed files with 996 additions and 991 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Date: 2023-06-23
Date: 2023-06-24
Package: CHNOSZ
Version: 2.0.0-12
Version: 2.0.0-13
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors@R: c(
person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),
Expand Down
14 changes: 7 additions & 7 deletions R/Berman.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,14 @@
# in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2.
# J. Petrol. 29, 445-522. https://doi.org/10.1093/petrology/29.2.445

Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE, calc.disorder=TRUE) {
Berman <- function(name, T = 298.15, P = 1, check.G = FALSE, calc.transition = TRUE, calc.disorder = TRUE) {
# Reference temperature and pressure
Pr <- 1
Tr <- 298.15
# Make T and P the same length
ncond <- max(length(T), length(P))
T <- rep(T, length.out=ncond)
P <- rep(P, length.out=ncond)
T <- rep(T, length.out = ncond)
P <- rep(P, length.out = ncond)

# Get parameters in the Berman equations
# Start with thermodynamic parameters provided with CHNOSZ
Expand All @@ -20,10 +20,10 @@ Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE,
userfile <- get("thermo", CHNOSZ)$opt$Berman
userfileexists <- FALSE
if(!is.na(userfile)) {
if(userfile!="") {
if(userfile != "") {
if(file.exists(userfile)) {
userfileexists <- TRUE
BDat_user <- read.csv(userfile, as.is=TRUE)
BDat_user <- read.csv(userfile, as.is = TRUE)
dat <- rbind(BDat_user, dat)
} else stop("the file named in thermo()$opt$Berman (", userfile, ") does not exist")
}
Expand All @@ -38,7 +38,7 @@ Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE,
if(missing(name)) return(dat) else {
# Which row has data for this mineral?
irow <- which(dat$name == name)
if(length(irow)==0) {
if(length(irow) == 0) {
if(userfileexists) stop("Data for ", name, " not available. Please add it to ", userfile)
if(!userfileexists) stop("Data for ", name, " not available. Please add it to your_data_file.csv and run thermo('opt$Berman' = 'path/to/your_data_file.csv')")
}
Expand All @@ -62,7 +62,7 @@ Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE,
Gdiff <- GfPrTr_calc - GfPrTr
#if(is.na(GfPrTr)) warning(paste0(name, ": GfPrTr(table) is NA"), call.=FALSE)
if(!is.na(GfPrTr)) if(abs(Gdiff) >= 1000) warning(paste0(name, ": GfPrTr(calc) - GfPrTr(table) is too big! == ",
round(GfPrTr_calc - GfPrTr), " J/mol"), call.=FALSE)
round(GfPrTr_calc - GfPrTr), " J/mol"), call. = FALSE)
# (the tabulated GfPrTr is unused below)
}

Expand Down
10 changes: 5 additions & 5 deletions R/DEW.R
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,8 @@ calculateDensity <- function(pressure, temperature, error = 0.01) {
calculateDensity
}
# 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))
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
sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
}
Expand Down Expand Up @@ -88,8 +88,8 @@ calculateGibbsOfWater <- function(pressure, temperature) {
GAtOneKb + integral
}
# 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))
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
sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
}
Expand Down Expand Up @@ -223,7 +223,7 @@ calculateOmega <- function(pressure, temperature, density, wref, Z) {
omega <- eta * (Z * Z / re - Z / (3.082 + g))
# 'If species is hydrogen, the species is neutral, or the pressure is above 6 kb,
# 'this equation is not necessary because omega is very close to wref.
if(Z==0) omega[] <- wref
if(Z == 0) omega[] <- wref
omega[pressure > 6000] <- wref
}

Expand Down
4 changes: 2 additions & 2 deletions R/IAPWS95.R
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Calculate properties of water using the IAPWS-95 formulation (Wagner and Pruss, 2002)

IAPWS95 <- function(property,T=298.15,rho=1000) {
IAPWS95 <- function(property,T = 298.15,rho = 1000) {
property <- tolower(property)
# Triple point
T.triple <- 273.16 # K
Expand Down Expand Up @@ -95,7 +95,7 @@ IAPWS95 <- function(property,T=298.15,rho=1000) {
t <- c(t,get(property[j])())
}
t <- data.frame(t)
if(j==1) ww <- t else ww <- cbind(ww,t)
if(j == 1) ww <- t else ww <- cbind(ww,t)
}
colnames(ww) <- property
return(ww)
Expand Down
27 changes: 14 additions & 13 deletions R/add.OBIGT.R
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,15 @@ mod.OBIGT <- function(..., zap = FALSE) {
if(is.list(args[[1]])) args <- args[[1]]
if(length(args) < 2) stop("please supply at least a species name and a property to update")
if(is.null(names(args))) stop("all arguments after the first should be named")
if(any(tail(nchar(names(args)), -1)==0)) stop("all arguments after the first should be named")
if(any(tail(nchar(names(args)), -1) == 0)) stop("all arguments after the first should be named")
# If the first argument is numeric, it's the species index
if(is.numeric(args[[1]][1])) {
ispecies <- args[[1]]
args <- args[-1]
speciesname <- info(ispecies, check.it = FALSE)$name
} else {
# If the name of the first argument is missing, assume it's the species name
if(names(args)[1]=="") names(args)[1] <- "name"
if(names(args)[1] == "") names(args)[1] <- "name"
speciesname <- args$name
# Search for this species, use check.protein = FALSE to avoid infinite loop when adding proteins
# and suppressMessages to not show messages about matches of this name to other states
Expand All @@ -33,18 +33,18 @@ mod.OBIGT <- function(..., zap = FALSE) {
species = args$name, check.protein = FALSE, SIMPLIFY = TRUE, USE.NAMES = FALSE))
}
# The column names of thermo()$OBIGT, split at the "."
cnames <- c(do.call(rbind, strsplit(colnames(thermo$OBIGT), ".", fixed=TRUE)), colnames(thermo$OBIGT))
cnames <- c(do.call(rbind, strsplit(colnames(thermo$OBIGT), ".", fixed = TRUE)), colnames(thermo$OBIGT))
# The columns we are updating
icol <- match(names(args), cnames)
if(any(is.na(icol))) stop(paste("properties not in thermo$OBIGT:", paste(names(args)[is.na(icol)], collapse=" ")) )
if(any(is.na(icol))) stop(paste("properties not in thermo$OBIGT:", paste(names(args)[is.na(icol)], collapse = " ")) )
# The column numbers for properties that matched after the split
icol[icol > 44] <- icol[icol > 44] - 44
icol[icol > 22] <- icol[icol > 22] - 22
# Which species are new and which are old
inew <- which(is.na(ispecies))
iold <- which(!is.na(ispecies))
# The arguments as data frame
args <- data.frame(args, stringsAsFactors=FALSE)
args <- data.frame(args, stringsAsFactors = FALSE)
if(length(inew) > 0) {
# The right number of blank rows of thermo()$OBIGT
newrows <- thermo$OBIGT[1:length(inew), ]
Expand All @@ -62,7 +62,7 @@ mod.OBIGT <- function(..., zap = FALSE) {
namodel <- is.na(newrows$model)
if(any(namodel)) newrows$model[namodel] <- ifelse(newrows$state[namodel] == "aq", "HKF", "CGL")
# Now check the formulas
e <- tryCatch(makeup(newrows$formula), error=function(e) e)
e <- tryCatch(makeup(newrows$formula), error = function(e) e)
if(inherits(e, "error")) {
warning("please supply a valid chemical formula as the species name or in the 'formula' argument")
# Transmit the error from makeup
Expand All @@ -76,7 +76,8 @@ mod.OBIGT <- function(..., zap = FALSE) {
ntotal <- nrow(thermo$OBIGT)
ispecies[inew] <- (ntotal-length(inew)+1):ntotal
# Inform user
message(paste("mod.OBIGT: added ", newrows$name, "(", newrows$state, ")", " with ", newrows$model, " model and energy units of ", newrows$E_units, sep="", collapse="\n"))
message(paste("mod.OBIGT: added ", newrows$name, "(", newrows$state, ")", " with ", newrows$model,
" model and energy units of ", newrows$E_units, sep = "", collapse = "\n"))
}
if(length(iold) > 0) {
# Loop over species
Expand Down Expand Up @@ -108,26 +109,26 @@ mod.OBIGT <- function(..., zap = FALSE) {
return(ispecies)
}

add.OBIGT <- function(file, species=NULL, force=TRUE) {
add.OBIGT <- function(file, species = NULL, force = TRUE) {
# Add/replace entries in thermo$OBIGT from values saved in a file
# Only replace if force==TRUE
# Only replace if force == TRUE
thermo <- get("thermo", CHNOSZ)
to1 <- thermo$OBIGT
id1 <- paste(to1$name,to1$state)
# We match system files with the file suffixes (.csv) removed
sysfiles <- dir(system.file("extdata/OBIGT/", package="CHNOSZ"))
sysfiles <- dir(system.file("extdata/OBIGT/", package = "CHNOSZ"))
sysnosuffix <- sapply(strsplit(sysfiles, "\\."), "[", 1)
isys <- match(file, sysnosuffix)
if(!is.na(isys)) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package="CHNOSZ")
if(!is.na(isys)) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package = "CHNOSZ")
# Read data from the file
to2 <- read.csv(file, as.is=TRUE)
to2 <- read.csv(file, as.is = TRUE)
Etxt <- paste(unique(to2$E_units), collapse = " and ")
# Load only selected species if requested
if(!is.null(species)) {
idat <- match(species, to2$name)
ina <- is.na(idat)
if(!any(ina)) to2 <- to2[idat, ]
else stop(paste("file", file, "doesn't have", paste(species[ina], collapse=", ")))
else stop(paste("file", file, "doesn't have", paste(species[ina], collapse = ", ")))
}
id2 <- paste(to2$name,to2$state)
# Check if the data table is compatible with thermo$OBIGT
Expand Down
2 changes: 1 addition & 1 deletion R/add.protein.R
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ add.protein <- function(aa, as.residue = FALSE) {
aasum <- function(aa, abundance = 1, average = FALSE, protein = NULL, organism = NULL) {
# Returns the sum of the amino acid counts in aa,
# multiplied by the abundances of the proteins
abundance <- rep(abundance, length.out=nrow(aa))
abundance <- rep(abundance, length.out = nrow(aa))
# Drop any NA rows or abundances
ina.aa <- is.na(aa$chains)
ina.ab <- is.na(abundance)
Expand Down
54 changes: 27 additions & 27 deletions R/affinity.R
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@
#source("util.data.R")
#source("species.R")

affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rhomin=FALSE,
return.buffer=FALSE, return.sout=FALSE, balance="PBB", iprotein=NULL, loga.protein=0, transect = NULL) {
affinity <- function(..., property = NULL, sout = NULL, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
return.buffer = FALSE, return.sout = FALSE, balance = "PBB", iprotein = NULL, loga.protein = 0, transect = NULL) {
# ...: variables over which to calculate
# property: what type of energy
# (G.basis,G.species,logact.basis,logK,logQ,A)
Expand All @@ -34,7 +34,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
args.orig <- list(...)
# We can only do anything with at least one argument
if(length(args.orig) > 0) {
if(identical(args.orig[[1]][1], list(fun="affinity"))) {
if(identical(args.orig[[1]][1], list(fun = "affinity"))) {
aargs <- args.orig[[1]]$args
# We can only update arguments given after the first argument
if(length(args.orig) > 1) {
Expand All @@ -49,7 +49,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho

# The argument list
args <- energy.args(args.orig, transect = transect)
args <- c(args, list(sout=sout, exceed.Ttr=exceed.Ttr, exceed.rhomin=exceed.rhomin))
args <- c(args, list(sout = sout, exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin))

# The species we're given
thermo <- get("thermo", CHNOSZ)
Expand Down Expand Up @@ -77,7 +77,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
# Add protein residues to the species list
resnames <- c("H2O",aminoacids(3))
# Residue activities set to zero; account for protein activities later
resprot <- paste(resnames,"RESIDUE",sep="_")
resprot <- paste(resnames,"RESIDUE",sep = "_")
species(resprot, 0)
thermo <- get("thermo", CHNOSZ)
ires <- match(resprot, thermo$species$name)
Expand All @@ -93,7 +93,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
buffer <- TRUE
message('affinity: loading buffer species')
if(!is.null(thermo$species)) is.species <- 1:nrow(thermo$species) else is.species <- numeric()
is.buffer <- buffer(logK=NULL)
is.buffer <- buffer(logK = NULL)
thermo <- get("thermo", CHNOSZ)
is.buff <- numeric()
for(i in 1:length(is.buffer)) is.buff <- c(is.buff,as.numeric(is.buffer[[i]]))
Expand All @@ -117,7 +117,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
logact.basis.new <- logact.basis <- do.call("energy",args)$a
ibasis.new <- numeric()
for(k in 1:length(buffers)) {
ibasis <- which(as.character(mybasis$logact)==buffers[k])
ibasis <- which(as.character(mybasis$logact) == buffers[k])
# Calculate the logKs from the affinities
logK <- a
for(i in 1:length(logK)) {
Expand All @@ -126,8 +126,8 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
logK[[i]] <- logK[[i]] - logact.basis.new[[j]] * thermo$species[i,j]
}
}
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)
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
is.only.buffer.new <- is.only.buffer[is.only.buffer %in% is.buffer[[k]]]
Expand All @@ -143,12 +143,12 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
#if(!identical(a[[i]],aa)) print(paste(i,j))
}
}
if(k==length(buffers) & return.buffer) {
if(k == length(buffers) & return.buffer) {
logact.basis.new <- lbn[[2]]
ibasis.new <- c(ibasis.new,lbn[[1]])
} else ibasis.new <- c(ibasis.new,ibasis)
}
species(is.only.buffer,delete=TRUE)
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
tb <- logact.basis.new[unique(ibasis.new)]
Expand All @@ -159,9 +159,9 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
for(i in 1:length(tb)) {
#tb[[i]] <- as.data.frame(tb[[i]])
if(nd > 0) rownames(tb[[i]]) <-
seq(args$lims[[1]][1],args$lims[[1]][2],length.out=args$lims[[1]][3])
seq(args$lims[[1]][1],args$lims[[1]][2],length.out = args$lims[[1]][3])
if(nd > 1) colnames(tb[[i]]) <-
seq(args$lims[[2]][1],args$lims[[2]][2],length.out=args$lims[[2]][3])
seq(args$lims[[2]][1],args$lims[[2]][2],length.out = args$lims[[2]][3])
}
}
}
Expand All @@ -171,7 +171,7 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
if(!is.null(iprotein)) {
# Fast protein calculations 20090331
# Function to calculate affinity of formation reactions from those of residues
loga.protein <- rep(loga.protein,length.out=length(iprotein))
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])
Expand All @@ -186,14 +186,14 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
# The current species list, containing the residues
resspecies <- thermo$species
# Now we can delete the residues from the species list
species(ires,delete=TRUE)
species(ires,delete = TRUE)
# State and protein names
state <- resspecies$state[1]
name <- paste(thermo$protein$protein[iprotein],thermo$protein$organism[iprotein],sep="_")
name <- paste(thermo$protein$protein[iprotein],thermo$protein$organism[iprotein],sep = "_")
# 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
protspecies <- cbind(protbasis,data.frame(ispecies=ispecies,logact=loga.protein,state=state,name=name))
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
Expand Down Expand Up @@ -231,8 +231,8 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
# they got converted to log_a(e-) at an unknown temperature
Eharg <- args.orig[[iEh]]
if(length(Eharg) > 3) Ehvals <- Eharg
else if(length(Eharg) == 3) Ehvals <- seq(Eharg[1], Eharg[2], length.out=Eharg[3])
else if(length(Eharg) == 2) Ehvals <- seq(Eharg[1], Eharg[2], length.out=256)
else if(length(Eharg) == 3) Ehvals <- seq(Eharg[1], Eharg[2], length.out = Eharg[3])
else if(length(Eharg) == 2) Ehvals <- seq(Eharg[1], Eharg[2], length.out = 256)
vals[[iEh]] <- Ehvals
}
# Get pe and pH
Expand All @@ -252,18 +252,18 @@ affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rho
names(vals) <- vars

# Content of return value depends on buffer request
if(return.buffer) return(c(tb, list(vars=vars, vals=vals)))
if(return.buffer) return(c(tb, list(vars = vars, vals = vals)))
# 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))
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
# (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,
basis=mybasis, species=myspecies, T=T, P=P, IS=args$IS, vars=vars, vals=vals, values=a)
else a <- list(fun="affinity", args=allargs, sout=sout, property=property,
basis=mybasis, species=myspecies, T=T, P=P, vars=vars, vals=vals, values=a)
if(buffer) a <- c(a, list(buffer=tb))
if(!is.na(iIS)) a <- list(fun = "affinity", args = allargs, sout = sout, property = property,
basis = mybasis, species = myspecies, T = T, P = P, IS = args$IS, vars = vars, vals = vals, values = a)
else a <- list(fun = "affinity", args = allargs, sout = sout, property = property,
basis = mybasis, species = myspecies, T = T, P = P, vars = vars, vals = vals, values = a)
if(buffer) a <- c(a, list(buffer = tb))
return(a)
}

Loading

0 comments on commit aa13309

Please sign in to comment.