Split original git repo by subdirectory.
Taha Ahmed
13 years ago
parent
5ee85d580e
commit
a6466ddc02
@ -1,7 +0,0 @@
|
||||
# AUTOLAB.R
|
||||
# Functions to read and manipulate data from AUTOLAB potentiostats/galvanostats
|
||||
# Taha Ahmed, June 2011
|
||||
|
||||
# CONTENTS
|
||||
# >>>> amp2df
|
||||
# >>>> ocp2df
|
||||
@ -1,89 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
#################### amp2df ######################
|
||||
##################################################
|
||||
amp2df <- function(datafilename, wearea = 1) {
|
||||
## Description:
|
||||
## Reads current-time data (from AUTOLAB potentiostat)
|
||||
## and returns a dataframe with the data and some
|
||||
## calculated quantities based on the data.
|
||||
## Usage:
|
||||
## amp2df(datafilename, wearea)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## wearea: (optional) area of working electrode (in square centimeters)
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ time : num
|
||||
## $ current : num
|
||||
## $ currentdensity : num
|
||||
## $ timediff : num
|
||||
## $ dIdt : num
|
||||
## $ didt : num
|
||||
## $ charge : num
|
||||
## $ chargedensity : num
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\d+\\.\\d+"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers _without_ a negative sign (hyphen),
|
||||
# followed by one or more digits before the decimal, a decimal point,
|
||||
# and an one or more digits after the decimal point.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(sampleid,
|
||||
stringsAsFactors = FALSE,
|
||||
matrix(scan(zz, what = numeric(), sep = ""),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "time", "current")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff, currentdensity = currentdensity)
|
||||
# Calculate time diff and current diff
|
||||
timediff <- c(ff$time[1], diff(ff$time))
|
||||
currentdiff <- c(ff$current[1], diff(ff$current))
|
||||
currentdensitydiff <- c(ff$currentdensity[1], diff(ff$currentdensity))
|
||||
# Calculate differential of current and current density
|
||||
dIdt <- currentdiff / timediff
|
||||
didt <- currentdensitydiff / timediff
|
||||
# Calculate charge and charge density
|
||||
charge <- cumsum(ff$current)
|
||||
chargedensity <- cumsum(ff$currentdensity)
|
||||
# Update ff dataframe
|
||||
ff <- cbind(ff,
|
||||
timediff = timediff,
|
||||
dIdt = dIdt,
|
||||
didt = didt,
|
||||
charge = charge,
|
||||
chargedensity = chargedensity)
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,61 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
#################### ocp2df ######################
|
||||
##################################################
|
||||
ocp2df <- function(datafilename) {
|
||||
## Description:
|
||||
## Reads voltage-time data (from AUTOLAB potentiostat)
|
||||
## and returns a dataframe with the data.
|
||||
## Usage:
|
||||
## ocp2df(datafilename)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ time : num
|
||||
## $ potential : num
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\d+\\.\\d+"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers _without_ a negative sign (hyphen),
|
||||
# followed by one or more digits before the decimal, a decimal point,
|
||||
# and an one or more digits after the decimal point.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(sampleid,
|
||||
stringsAsFactors = FALSE,
|
||||
matrix(scan(zz, what = numeric(), sep = ""),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "time", "potential")
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,12 +0,0 @@
|
||||
# CHI.R
|
||||
# Functions to read and manipulate data from the CHI760 potentiostat/galvanostat
|
||||
# Taha Ahmed, Jan 2011 - Feb 2011
|
||||
|
||||
# CONTENTS
|
||||
# >>>> mps2df
|
||||
# >>>> ocp2df
|
||||
# >>>> chronocm2df
|
||||
# >>>> chronoamp2df
|
||||
# >>>> amperometry2df
|
||||
# >>>> cv2df
|
||||
# >>>> lsv2df
|
||||
@ -1,123 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
############### amperometry2df ###################
|
||||
##################################################
|
||||
amperometry2df <- function(datafilename, wearea = 1) {
|
||||
## Description:
|
||||
## Reads current-time data (from CHI 760 potentiostat)
|
||||
## and returns a dataframe with the data,
|
||||
## the data attributes (experimental conditions),
|
||||
## and some calculated parameters (charge, didt, etc.)
|
||||
## Usage:
|
||||
## amperometry2df(datafilename, wearea)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## wearea: (optional) area of working electrode (in square centimeter)
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr
|
||||
## $ time : num
|
||||
## $ current : num
|
||||
## $ currentdensity : num
|
||||
## $ timediff : num
|
||||
## $ dIdt : num
|
||||
## $ didt : num
|
||||
## $ charge : num
|
||||
## $ chargedensity : num
|
||||
## $ InitE : num
|
||||
## $ SampleInterval : num
|
||||
## $ RunTime : num
|
||||
## $ QuietTime : num
|
||||
## $ Sensitivity : num
|
||||
## Note:
|
||||
## The CH Instruments 760 potentiostat records all data
|
||||
## using standard SI units, therefore this function
|
||||
## assumes all potential values to be in volts,
|
||||
## currents to be in amperes, charges in Coulombs,
|
||||
## time in seconds, and so on.
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+[e,]"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by either the letter 'e' or a comma.
|
||||
# Note that backslashes are escaped due to the way R handles strings.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(stringsAsFactors = FALSE,
|
||||
sampleid, matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "time", "current")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff, currentdensity = currentdensity)
|
||||
# Calculate time and current diffs
|
||||
timediff <- c(ff$time[1], diff(ff$time))
|
||||
currentdiff <- c(ff$current[1], diff(ff$current))
|
||||
currentdensitydiff <- c(ff$currentdensity[1], diff(ff$currentdensity))
|
||||
# Calculate differential of current and current density
|
||||
dIdt <- currentdiff / timediff
|
||||
didt <- currentdensitydiff / timediff
|
||||
# Calculate charge and charge density
|
||||
charge <- cumsum(ff$current)
|
||||
chargedensity <- cumsum(ff$currentdensity)
|
||||
# Update ff dataframe
|
||||
ff <- cbind(ff,
|
||||
timediff = timediff,
|
||||
dIdt = dIdt,
|
||||
didt = didt,
|
||||
charge = charge,
|
||||
chargedensity = chargedensity)
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# InitE (volt)
|
||||
position.InitE <- regexpr("^Init\\sE\\s\\(V\\)", chifile)
|
||||
InitE <- as.numeric(strsplit(chifile[which(position.InitE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$InitE <- InitE
|
||||
# SampleInterval (volt)
|
||||
position.SampleInterval <- regexpr("^Sample\\sInterval\\s\\(s\\)", chifile)
|
||||
SampleInterval <- as.numeric(strsplit(chifile[which(position.SampleInterval == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$SampleInterval <- SampleInterval
|
||||
# Run time (seconds)
|
||||
position.RunTime <- regexpr("^Run\\sTime\\s\\(sec\\)", chifile)
|
||||
RunTime <- as.numeric(strsplit(chifile[which(position.RunTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$RunTime <- RunTime
|
||||
# Quiet time (seconds)
|
||||
position.QuietTime <- regexpr("^Quiet\\sTime\\s\\(sec\\)", chifile)
|
||||
QuietTime <- as.numeric(strsplit(chifile[which(position.QuietTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$QuietTime <- QuietTime
|
||||
# Sensitivity (ampere per volt)
|
||||
position.Sensitivity <- regexpr("^Sensitivity\\s\\(A/V\\)", chifile)
|
||||
Sensitivity <- as.numeric(strsplit(chifile[which(position.Sensitivity == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Sensitivity <- Sensitivity
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,122 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
################# chronoamp2df ###################
|
||||
##################################################
|
||||
chronoamp2df <- function(datafilename, wearea = 1) {
|
||||
## Description:
|
||||
## Reads current-time data (from CHI 760 potentiostat)
|
||||
## and returns a dataframe with the data and the
|
||||
## data attributes (experimental conditions).
|
||||
## Usage:
|
||||
## chronoamp2df(datafilename, wearea)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## wearea: (optional) area of working electrode (in square centimeters)
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr
|
||||
## $ step : num
|
||||
## $ time : num
|
||||
## $ current : num
|
||||
## $ currentdensity : num
|
||||
## $ InitE : num
|
||||
## $ HighE : num
|
||||
## $ LowE : num
|
||||
## $ InitPN : chr
|
||||
## $ Step : num
|
||||
## $ Pulse width : num
|
||||
## $ Sample interval : num
|
||||
## $ Quiet Time : num
|
||||
## $ Sensitivity : num
|
||||
## Note:
|
||||
## The CH Instruments 760 potentiostat records all data
|
||||
## using standard SI units, therefore this function
|
||||
## assumes all potential values to be in volts,
|
||||
## currents to be in amperes, charges in Coulombs,
|
||||
## time in seconds, and so on.
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+[e,]"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by either the letter 'e' or a comma.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(sampleid, step.current = s,
|
||||
matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "step.current", "time", "current")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff, currentdensity = currentdensity)
|
||||
#
|
||||
|
||||
### Collect attributes of this experiment
|
||||
# InitE (volt)
|
||||
position.InitE <- regexpr("^Init\\sE\\s\\(V\\)", chifile)
|
||||
InitE <- as.numeric(strsplit(chifile[which(position.InitE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$InitE <- InitE
|
||||
# HighE (volt)
|
||||
position.HighE <- regexpr("^High\\sE\\s\\(V\\)", chifile)
|
||||
HighE <- as.numeric(strsplit(chifile[which(position.HighE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$HighE <- HighE
|
||||
# LowE (volt)
|
||||
position.LowE <- regexpr("^Low\\sE\\s\\(V\\)", chifile)
|
||||
LowE <- as.numeric(strsplit(chifile[which(position.LowE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$LowE <- LowE
|
||||
# InitPN
|
||||
position.InitPN <- regexpr("^Init\\sP/N", chifile)
|
||||
InitPN <- strsplit(chifile[which(position.InitPN == 1)], "\\s=\\s")[[1]][2]
|
||||
ff$InitPN <- InitPN
|
||||
# Steps (total number of steps)
|
||||
position.Steps <- regexpr("^Step\\s", chifile)
|
||||
Steps <- as.numeric(strsplit(chifile[which(position.Steps == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Steps <- Steps
|
||||
# Pulse width (s)
|
||||
position.PulseWidth <- regexpr("^Pulse\\sWidth\\s\\(sec\\)", chifile)
|
||||
PulseWidth <- as.numeric(strsplit(chifile[which(position.PulseWidth == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$PulseWidth <- PulseWidth
|
||||
# Sample interval (s)
|
||||
position.SampleInterval <- regexpr("^Sample\\sInterval\\s\\(s\\)", chifile)
|
||||
SampleInterval <- as.numeric(strsplit(chifile[which(position.SampleInterval == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$SampleInterval <- SampleInterval
|
||||
# Quiet Time (s)
|
||||
position.QuietTime <- regexpr("^Quiet\\sTime\\s\\(sec\\)", chifile)
|
||||
QuietTime <- as.numeric(strsplit(chifile[which(position.QuietTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$QuietTime <- QuietTime
|
||||
# Sensitivity (A/V)
|
||||
position.Sensitivity <- regexpr("^Sensitivity\\s\\(A/V\\)", chifile)
|
||||
Sensitivity <- as.numeric(strsplit(chifile[which(position.Sensitivity == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Sensitivity <- Sensitivity
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,78 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
################# chronocm2df ####################
|
||||
##################################################
|
||||
chronocm2df <- function(datafilename) {
|
||||
# Function description: chronocoulometry data
|
||||
# CH Instruments potentiostat records all data using standard SI units,
|
||||
# so all potential values are in volts, currents are in amperes,
|
||||
# charges in Coulombs, time in seconds, etc.
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+[e,]"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by either the letter 'e' or a comma.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(sampleid, step = factor(s),
|
||||
matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "step", "time", "charge")
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# These attributes are specific for each kind of experiment,
|
||||
# be careful when adapting to other electrochemical data
|
||||
rgxp.attr <- c("^Init\\sE\\s\\(V\\)",
|
||||
"^Final\\sE\\s\\(V\\)",
|
||||
"^Step\\s",
|
||||
"^Pulse\\sWidth\\s\\(sec\\)",
|
||||
"^Sample\\sInterval\\s\\(s\\)",
|
||||
"^Quiet\\sTime\\s\\(sec\\)",
|
||||
"^Sensitivity\\s\\(A/V\\)")
|
||||
names.attr <- c("InitE",
|
||||
"FinalE",
|
||||
"Steps",
|
||||
"PulseWidth",
|
||||
"SamplingInterval",
|
||||
"QuietTime",
|
||||
"Sensitivity")
|
||||
for (n in 1:length(rgxp.attr)) {
|
||||
attrow.idx <- regexpr(rgxp.attr[n], chifile)
|
||||
attrow.len <- attr(attrow.idx, "match.length")
|
||||
attr(attrow.idx, "match.length") <- NULL
|
||||
attr(ff, names.attr[n]) <- strsplit(chifile[which(attrow.idx == 1)],
|
||||
"\\s=\\s")[[1]][2]
|
||||
}
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,94 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
#################### cv2df #######################
|
||||
##################################################
|
||||
cv2df <- function(datafilename, wearea = 1) {
|
||||
# Function description:
|
||||
# CH Instruments potentiostat records all data using standard SI units,
|
||||
# so all potential values are in volts, currents are in amperes,
|
||||
# charges in Coulombs, time in seconds, etc.
|
||||
#
|
||||
cvfile <- file(datafilename, "r")
|
||||
chifile <- readLines(cvfile, n = -1) #read all lines of input file
|
||||
close(cvfile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+,"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by a comma.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(sampleid, cycle = as.integer(ceiling(s/2)), segment = s,
|
||||
matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 3, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
# Column names after initial assignment
|
||||
names(ff) <- c("sampleid", "cycle", "segment", "potential", "current", "charge")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff[, 1:5], currentdensity = currentdensity, charge = ff[, 6])
|
||||
## Collect attributes of this experiment
|
||||
# InitE (volt)
|
||||
position.InitE <- regexpr("^Init\\sE\\s\\(V\\)", chifile)
|
||||
InitE <- as.numeric(strsplit(chifile[which(position.InitE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$InitE <- InitE
|
||||
# HighE (volt)
|
||||
position.HighE <- regexpr("^High\\sE\\s\\(V\\)", chifile)
|
||||
HighE <- as.numeric(strsplit(chifile[which(position.HighE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$HighE <- HighE
|
||||
# LowE (volt)
|
||||
position.LowE <- regexpr("^Low\\sE\\s\\(V\\)", chifile)
|
||||
LowE <- as.numeric(strsplit(chifile[which(position.LowE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$LowE <- LowE
|
||||
# InitPN (positive or negative)
|
||||
position.InitPN <- regexpr("^Init\\sP/N", chifile)
|
||||
InitPN <- strsplit(chifile[which(position.InitPN == 1)], "\\s=\\s")[[1]][2]
|
||||
ff$InitPN <- InitPN
|
||||
# ScanRate (volt per second)
|
||||
position.ScanRate <- regexpr("^Scan\\sRate\\s\\(V/s\\)", chifile)
|
||||
ScanRate <- as.numeric(strsplit(chifile[which(position.ScanRate == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$ScanRate <- ScanRate
|
||||
# Segments, number of
|
||||
position.Segments <- regexpr("^Segment\\s=", chifile)
|
||||
Segments <- as.numeric(strsplit(chifile[which(position.Segments == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Segments <- Segments
|
||||
# SampleInterval (volt)
|
||||
position.SampleInterval <- regexpr("^Sample\\sInterval\\s\\(V\\)", chifile)
|
||||
SampleInterval <- as.numeric(strsplit(chifile[which(position.SampleInterval == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$SampleInterval <- SampleInterval
|
||||
# Quiet time (seconds)
|
||||
position.QuietTime <- regexpr("^Quiet\\sTime\\s\\(sec\\)", chifile)
|
||||
QuietTime <- as.numeric(strsplit(chifile[which(position.QuietTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$QuietTime <- QuietTime
|
||||
# Sensitivity (ampere per volt)
|
||||
position.Sensitivity <- regexpr("^Sensitivity\\s\\(A/V\\)", chifile)
|
||||
Sensitivity <- as.numeric(strsplit(chifile[which(position.Sensitivity == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Sensitivity <- Sensitivity
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,114 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
################### lsv2df #######################
|
||||
##################################################
|
||||
lsv2df <- function(datafilename, wearea = 1) {
|
||||
## Description:
|
||||
## Reads LSV datafiles from CHI 760 potentiostat
|
||||
## (potential, current, and charge)
|
||||
## and returns a dataframe with the data,
|
||||
## the data attributes (experimental conditions),
|
||||
## and calculated current density and charge density.
|
||||
## Usage:
|
||||
## lsv2df(datafilename, wearea)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## wearea: (optional) area of working electrode (in square centimeter)
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr (id)
|
||||
## $ segment : num (id)
|
||||
## $ potential : num (measure)
|
||||
## $ current : num (measure)
|
||||
## $ charge : num (measure)
|
||||
## $ currentdensity : num (measure)
|
||||
## $ chargedensity : num (measure)
|
||||
## $ InitE : num (id)
|
||||
## $ FinalE : num (id)
|
||||
## $ ScanRate : num (id)
|
||||
## $ SampleInterval : num (id)
|
||||
## $ QuietTime : num (id)
|
||||
## $ Sensitivity : num (id)
|
||||
## Note:
|
||||
## The CH Instruments 760 potentiostat records all data
|
||||
## using standard SI units, therefore this function
|
||||
## assumes all potential values to be in volts,
|
||||
## currents to be in amperes, charges in Coulombs,
|
||||
## time in seconds, and so on.
|
||||
#
|
||||
lsvfile <- file(datafilename, "r")
|
||||
chifile <- readLines(lsvfile, n = -1) #read all lines of input file
|
||||
close(lsvfile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+,"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by a comma.
|
||||
# Note that backslashes are escaped.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(stringsAsFactors = FALSE,
|
||||
sampleid, segment = s,
|
||||
matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 3, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "segment", "potential", "current", "charge")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff, currentdensity = currentdensity)
|
||||
# Calculate charge density
|
||||
chargedensity <- ff$charge / wearea
|
||||
ff <- cbind(ff, chargedensity = chargedensity)
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# InitE (volt)
|
||||
position.InitE <- regexpr("^Init\\sE\\s\\(V\\)", chifile)
|
||||
InitE <- as.numeric(strsplit(chifile[which(position.InitE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$InitE <- InitE
|
||||
# FinalE (volt)
|
||||
position.FinalE <- regexpr("^Final\\sE\\s\\(V\\)", chifile)
|
||||
FinalE <- as.numeric(strsplit(chifile[which(position.FinalE == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$FinalE <- FinalE
|
||||
# ScanRate (volt per second)
|
||||
position.ScanRate <- regexpr("^Scan\\sRate\\s\\(V/s\\)", chifile)
|
||||
ScanRate <- as.numeric(strsplit(chifile[which(position.ScanRate == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$ScanRate <- ScanRate
|
||||
# SampleInterval (volt)
|
||||
position.SampleInterval <- regexpr("^Sample\\sInterval\\s\\(V\\)", chifile)
|
||||
SampleInterval <- as.numeric(strsplit(chifile[which(position.SampleInterval == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$SampleInterval <- SampleInterval
|
||||
# Quiet time (seconds)
|
||||
position.QuietTime <- regexpr("^Quiet\\sTime\\s\\(sec\\)", chifile)
|
||||
QuietTime <- as.numeric(strsplit(chifile[which(position.QuietTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$QuietTime <- QuietTime
|
||||
# Sensitivity (ampere per volt)
|
||||
position.Sensitivity <- regexpr("^Sensitivity\\s\\(A/V\\)", chifile)
|
||||
Sensitivity <- as.numeric(strsplit(chifile[which(position.Sensitivity == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Sensitivity <- Sensitivity
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,147 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
################### mps2df #######################
|
||||
##################################################
|
||||
mps2df <- function(datafilename, wearea = 1) {
|
||||
## Description:
|
||||
## Reads time vs current data (from CHI 760 potentiostat)
|
||||
## and returns a dataframe with the data and
|
||||
## the data attributes (experimental conditions).
|
||||
## Usage:
|
||||
## mps2df(datafilename, wearea)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## wearea: working electrode area in square centimeters (optional)
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr
|
||||
## $ time : num
|
||||
## $ current : num
|
||||
## $ currentdensity : num
|
||||
## $ timediff : num
|
||||
## $ dIdt : num
|
||||
## $ didt : num
|
||||
## $ charge : num
|
||||
## $ chargedensity : num
|
||||
## $ PotentialSteps : num
|
||||
## $ TimeSteps : num
|
||||
## $ Cycle : num
|
||||
## $ SampleInterval : num
|
||||
## $ QuietTime : num
|
||||
## $ Sensitivity : num
|
||||
## Note:
|
||||
## The CH Instruments 760 potentiostat records all data
|
||||
## using standard SI units, therefore this function
|
||||
## assumes all potential values to be in volts,
|
||||
## currents to be in amperes, charges in Coulombs,
|
||||
## time in seconds, and so on.
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d+\\.\\d+[e,]"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by either the letter 'e' or a comma.
|
||||
# Note that backslashes are escaped due to the way R handles strings.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(stringsAsFactors = FALSE,
|
||||
sampleid, matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "time", "current")
|
||||
# Calculate current density
|
||||
currentdensity <- ff$current / wearea
|
||||
ff <- cbind(ff, currentdensity = currentdensity)
|
||||
# Calculate time diff and current diff
|
||||
timediff <- c(ff$time[1], diff(ff$time))
|
||||
currentdiff <- c(ff$current[1], diff(ff$current))
|
||||
currentdensitydiff <- c(ff$currentdensity[1], diff(ff$currentdensity))
|
||||
# Calculate differential of current and current density
|
||||
dIdt <- currentdiff / timediff
|
||||
didt <- currentdensitydiff / timediff
|
||||
# Calculate charge and charge density
|
||||
charge <- cumsum(ff$current)
|
||||
chargedensity <- cumsum(ff$currentdensity)
|
||||
# Update ff dataframe
|
||||
ff <- cbind(ff,
|
||||
timediff = timediff,
|
||||
dIdt = dIdt,
|
||||
didt = didt,
|
||||
charge = charge,
|
||||
chargedensity = chargedensity)
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# Potential steps
|
||||
PotentialSteps <- ""
|
||||
positions.PotentialSteps <- regexpr("^E\\d+\\s\\(V\\)", chifile)
|
||||
PotentialSteps <- paste("\\num{", strsplit(chifile[which(positions.PotentialSteps == 1)],
|
||||
"\\s=\\s")[[1]][2], "}", sep = "")
|
||||
if (length(which(positions.PotentialSteps == 1)) > 1) {
|
||||
for (i in 2:length(which(positions.PotentialSteps == 1))) {
|
||||
PotentialSteps <-
|
||||
paste(PotentialSteps,
|
||||
paste("\\num{", strsplit(chifile[which(positions.PotentialSteps == 1)],
|
||||
"\\s=\\s")[[i]][2], "}", sep = ""), sep = ", ")
|
||||
}
|
||||
}
|
||||
ff$PotentialSteps <- PotentialSteps
|
||||
# Time steps
|
||||
TimeSteps <- ""
|
||||
positions.TimeSteps <- regexpr("^T\\d+\\s\\(s\\)", chifile)
|
||||
TimeSteps <- paste("\\num{", strsplit(chifile[which(positions.TimeSteps == 1)],
|
||||
"\\s=\\s")[[1]][2], "}", sep = "")
|
||||
if (length(which(positions.TimeSteps == 1)) > 1) {
|
||||
for (i in 2:length(which(positions.TimeSteps == 1))) {
|
||||
TimeSteps <-
|
||||
paste(TimeSteps,
|
||||
paste("\\num{", strsplit(chifile[which(positions.TimeSteps == 1)],
|
||||
"\\s=\\s")[[i]][2], "}", sep = ""), sep = ", ")
|
||||
}
|
||||
}
|
||||
ff$TimeSteps <- TimeSteps
|
||||
# Cycles
|
||||
position.Cycles <- regexpr("^Cycle", chifile)
|
||||
Cycles <- as.numeric(strsplit(chifile[which(position.Cycles == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Cycles <- Cycles
|
||||
# Sample interval (s)
|
||||
position.SampleInterval <- regexpr("^Sample\\sInterval\\s\\(s\\)", chifile)
|
||||
SampleInterval <- as.numeric(strsplit(chifile[which(position.SampleInterval == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$SampleInterval <- SampleInterval
|
||||
# Quiet Time (s)
|
||||
position.QuietTime <- regexpr("^Quiet\\sTime\\s\\(sec\\)", chifile)
|
||||
QuietTime <- as.numeric(strsplit(chifile[which(position.QuietTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$QuietTime <- QuietTime
|
||||
# Sensitivity (A/V)
|
||||
position.Sensitivity <- regexpr("^Sensitivity\\s\\(A/V\\)", chifile)
|
||||
Sensitivity <- as.numeric(strsplit(chifile[which(position.Sensitivity == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$Sensitivity <- Sensitivity
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,85 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
##################################################
|
||||
################### ocp2df #######################
|
||||
##################################################
|
||||
ocp2df <- function(datafilename) {
|
||||
## Description:
|
||||
## Reads time vs potential data (from CHI 760 potentiostat)
|
||||
## and returns a dataframe with the data and
|
||||
## the data attributes (experimental conditions).
|
||||
## Usage:
|
||||
## ocp2df(datafilename)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr
|
||||
## $ time : num
|
||||
## $ current : num
|
||||
## $ currentdensity : num
|
||||
## $ timediff : num
|
||||
## $ dIdt : num
|
||||
## $ didt : num
|
||||
## $ charge : num
|
||||
## $ chargedensity : num
|
||||
## $ InitE : num
|
||||
## $ SampleInterval : num
|
||||
## $ RunTime : num
|
||||
## $ QuietTime : num
|
||||
## $ Sensitivity : num
|
||||
## Note:
|
||||
## The CH Instruments 760 potentiostat records all data
|
||||
## using standard SI units, therefore this function
|
||||
## assumes all potential values to be in volts,
|
||||
## currents to be in amperes, charges in Coulombs,
|
||||
## time in seconds, and so on.
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
rgxp.number <- "^\\-?\\d\\.\\d+[e,]"
|
||||
# regexp that matches a decimal number at the beginning of the line.
|
||||
# Matches numbers with or without a negative sign (hyphen),
|
||||
# followed by one digit before the decimal, a decimal point,
|
||||
# and an arbitrary number of digits after the decimal point,
|
||||
# immediately followed by either the letter 'e' or a comma.
|
||||
# Note that backslashes are escaped due to the way R handles strings.
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.number, chifile)
|
||||
# Save the match length attribute to another variable,
|
||||
numrow.len <- attr(numrow.idx, "match.length")
|
||||
# then scrap the attribute of the original variable.
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start indices of data ranges
|
||||
starts <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1) + 1
|
||||
# End indices, except for the last
|
||||
ends <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1)
|
||||
# Fix the last index of end indices
|
||||
ends <- c(ends, length(numrow.idx))
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(chifile[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff,
|
||||
data.frame(stringsAsFactors = FALSE,
|
||||
sampleid, matrix(scan(zz, what = numeric(), sep = ","),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
}
|
||||
names(ff) <- c("sampleid", "time", "potential")
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# RunTime (sec)
|
||||
position.RunTime <- regexpr("^Run\\sTime\\s\\(sec\\)", chifile)
|
||||
RunTime <- as.numeric(strsplit(chifile[which(position.RunTime == 1)], "\\s=\\s")[[1]][2])
|
||||
ff$RunTime <- RunTime
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,8 +0,0 @@
|
||||
# INCA.R
|
||||
# Functions to read and manipulate data from the Oxford INCA EDS
|
||||
# Taha Ahmed, April 2011
|
||||
|
||||
# CONTENTS
|
||||
source("/home/taha/chepec/chetex/common/R/common.R")
|
||||
# >>>> eds2df
|
||||
# >>>> edspk
|
||||
@ -1,88 +0,0 @@
|
||||
## This function should be renamed to EMSA2df(). The file it reads is a so-called EMSA spectral data file.
|
||||
# /TA, 111110
|
||||
|
||||
eds2df <- function(edstxtfile) {
|
||||
## Description:
|
||||
## Reads EDS textfile from INCA EDS.
|
||||
## Stores data in data frame.
|
||||
## Usage:
|
||||
## eds2df(edstxtfile)
|
||||
## Arguments:
|
||||
## edstxtfile: character string, the full filename
|
||||
## (with path) to one txt file.
|
||||
## Value:
|
||||
## A dataframe
|
||||
#
|
||||
incatxt <- file(edstxtfile, "r")
|
||||
edsfile <- readLines(incatxt, n = -1) #read all lines of input file
|
||||
close(incatxt)
|
||||
#
|
||||
sampleid <- ProvideSampleId(edstxtfile)
|
||||
#
|
||||
rgxp.comment <- "^\\#"
|
||||
#
|
||||
numrow.idx <- regexpr(rgxp.comment, edsfile)
|
||||
# scrap the match.length attribute
|
||||
attr(numrow.idx, "match.length") <- NULL
|
||||
#
|
||||
i <- seq(1, length(numrow.idx) - 1, 1)
|
||||
j <- seq(2, length(numrow.idx), 1)
|
||||
# Start index of data range
|
||||
start.idx <- which(numrow.idx[i] == 1 & numrow.idx[j] != 1) + 1
|
||||
# End index of the data range
|
||||
end.idx <- which(numrow.idx[i] != 1 & numrow.idx[j] == 1)
|
||||
#
|
||||
zz <- textConnection(edsfile[start.idx:end.idx], "r")
|
||||
#
|
||||
ff <- data.frame()
|
||||
ff <- data.frame(stringsAsFactors = FALSE,
|
||||
sampleid = sampleid,
|
||||
matrix(scan(zz, what = numeric(), sep = ","), ncol = 2, byrow = T))
|
||||
close(zz)
|
||||
names(ff) <- c("sampleid", "energy", "counts")
|
||||
#
|
||||
### Collect attributes of this experiment
|
||||
# XUnit
|
||||
position.XUnit <- regexpr("^\\#XUNITS", edsfile)
|
||||
XUnit <- as.character(strsplit(edsfile[which(position.XUnit == 1)], ":\\s")[[1]][2])
|
||||
ff$XUnit <- XUnit
|
||||
# YUnit
|
||||
position.YUnit <- regexpr("^\\#YUNITS", edsfile)
|
||||
YUnit <- as.character(strsplit(edsfile[which(position.YUnit == 1)], ":\\s")[[1]][2])
|
||||
ff$YUnit <- YUnit
|
||||
# Date
|
||||
position.Date <- regexpr("^\\#DATE", edsfile)
|
||||
Date <- strsplit(edsfile[which(position.Date == 1)], ":\\s")[[1]][2]
|
||||
ff$Date <- Date
|
||||
# Time
|
||||
position.Time <- regexpr("^\\#TIME", edsfile)
|
||||
Time <- strsplit(edsfile[which(position.Time == 1)], ":\\s")[[1]][2]
|
||||
ff$Time <- Time
|
||||
# XPerChannel
|
||||
position.XPerChannel <- regexpr("^\\#XPERCHAN", edsfile)
|
||||
XPerChannel <- as.numeric(strsplit(edsfile[which(position.XPerChannel == 1)], ":\\s")[[1]][2])
|
||||
ff$XPerChannel <- XPerChannel
|
||||
# Offset
|
||||
position.Offset <- regexpr("^\\#OFFSET", edsfile)
|
||||
Offset <- as.numeric(strsplit(edsfile[which(position.Offset == 1)], ":\\s")[[1]][2])
|
||||
ff$Offset <- Offset
|
||||
# ChOffset
|
||||
position.ChOffset <- regexpr("^\\#CHOFFSET", edsfile)
|
||||
ChOffset <- as.numeric(strsplit(edsfile[which(position.ChOffset == 1)], ":\\s")[[1]][2])
|
||||
ff$ChOffset <- ChOffset
|
||||
# LiveTime
|
||||
position.LiveTime <- regexpr("^\\#LIVETIME", edsfile)
|
||||
LiveTime <- as.numeric(strsplit(edsfile[which(position.LiveTime == 1)], ":\\s")[[1]][2])
|
||||
ff$LiveTime <- LiveTime
|
||||
# DeadTime is calculated from: REALTIME - LIVETIME
|
||||
position.RealTime <- regexpr("^\\#REALTIME", edsfile)
|
||||
RealTime <- as.numeric(strsplit(edsfile[which(position.RealTime == 1)], ":\\s")[[1]][2])
|
||||
DeadTime <- RealTime - LiveTime
|
||||
ff$DeadTime <- DeadTime
|
||||
# BeamEnergy
|
||||
position.BeamEnergy <- regexpr("^\\#BEAMKV", edsfile)
|
||||
BeamEnergy <- as.numeric(strsplit(edsfile[which(position.BeamEnergy == 1)], ":\\s")[[1]][2])
|
||||
ff$BeamEnergy <- BeamEnergy
|
||||
#
|
||||
return(ff)
|
||||
}
|
||||
@ -1,63 +0,0 @@
|
||||
edsWrapper <-
|
||||
function(data.exp, run, override = FALSE,
|
||||
kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1, maxwdth=0.20) {
|
||||
|
||||
print("... Started edsWrapper")
|
||||
|
||||
# check if edspk has already completed successfully for the current job
|
||||
current.dirname <- getwd()
|
||||
print(current.dirname)
|
||||
current.filename <- "eds-peak-data.rda"
|
||||
edsdatafile <- paste(current.dirname, current.filename, sep = "/")
|
||||
|
||||
|
||||
if (file.exists(edsdatafile) && !override) {
|
||||
print("... Started if-clause 1")
|
||||
|
||||
# File already exists
|
||||
# return the data using load() or data()
|
||||
|
||||
load(file = edsdatafile)
|
||||
|
||||
if (run > length(edsres)) {
|
||||
|
||||
print("... Started if-clause 1:1")
|
||||
|
||||
# then it does not really exist
|
||||
edsres[[run]] <- edspk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
save(edsres, file = edsdatafile)
|
||||
|
||||
print("... Ended if-clause 1:1")
|
||||
}
|
||||
|
||||
print("... Ended if-clause 1")
|
||||
|
||||
return(edsres)
|
||||
} else {
|
||||
|
||||
print("... Started else-clause 1")
|
||||
|
||||
if (!exists("edsres")) {
|
||||
edsres <- list()
|
||||
print("... edsres list created")
|
||||
}
|
||||
|
||||
# Need to call edspk() and save its results to file as above
|
||||
edsres[[run]] <- edspk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
save(edsres, file = edsdatafile)
|
||||
|
||||
print("... Ended else-clause 1")
|
||||
|
||||
return(edsres)
|
||||
}
|
||||
}
|
||||
@ -1,92 +0,0 @@
|
||||
edspk <-
|
||||
function(eds.exp, kerpk = 1, fitmaxiter = 50, gam = 1.0, scl.factor = 0.1, maxwdth=0.20) {
|
||||
|
||||
eds.base <- baselinefit(eds.exp, tau=2.0, gam=gam, scl.factor=scl.factor, maxwdth=maxwdth)
|
||||
|
||||
# This loop deals with the output from baselinefit()
|
||||
# It makes a "melted" dataframe in long form for each
|
||||
# separated peak for some baseline parameters
|
||||
eds.pks <- data.frame()
|
||||
eds.pks.basl <- data.frame()
|
||||
eds.pks.pmg <- data.frame()
|
||||
eds.pks.spl <- data.frame()
|
||||
peaks <- 1:length(eds.base$npks)
|
||||
for (s in peaks) {
|
||||
# recorded data in long form by separated peak
|
||||
eds.pks <- rbind(eds.pks, # column names assigned after loop
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
eds.exp[eds.base$indlsep[s]:eds.base$indrsep[s], ]))
|
||||
# the calculated baseline in long form by separated peak
|
||||
eds.pks.basl <- rbind(eds.pks.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = eds.exp[eds.base$indlsep[s]:eds.base$indrsep[s]],
|
||||
y = eds.base$baseline$basisl[eds.base$indlsep[s]:eds.base$indrsep[s]]))
|
||||
# the taut string estimation in long form by separated peak
|
||||
eds.pks.pmg <- rbind(eds.pks.pmg,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = eds.exp[eds.base$indlsep[s]:eds.base$indrsep[s]],
|
||||
y = eds.base$pmg$fn[eds.base$indlsep[s]:eds.base$indrsep[s]]))
|
||||
# the weighted smoothed spline in long form by separated peak
|
||||
eds.pks.spl <- rbind(eds.pks.spl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = eds.exp[eds.base$indlsep[s]:eds.base$indrsep[s]],
|
||||
y = eds.base$spl$reg[eds.base$indlsep[s]:eds.base$indrsep[s]]))
|
||||
}
|
||||
# Column names assigned to d.pks
|
||||
names(eds.pks) <- c("peak", "kernel", "x", "y")
|
||||
|
||||
|
||||
# This loop calls pkdecompint() on each peak separately
|
||||
# It makes a "melted" dataframe in long form for:
|
||||
eds.fit <- list() # holds pkdecompint output
|
||||
eds.fit.fitpk <- data.frame() # contains fitting curves
|
||||
eds.fit.parpk <- data.frame() # physical parameters by peak and kernel
|
||||
eds.nobasl <- data.frame() # data with baseline removed
|
||||
peaks <- 1:length(eds.base$npks)
|
||||
for (s in peaks) {
|
||||
######## PKDECOMPINT ########
|
||||
if (length(kerpk) > 1) {
|
||||
# set number of kernels per peak manually
|
||||
eds.fit[[s]] <- pkdecompint(eds.base, intnum = s,
|
||||
k = kerpk[s], maxiter = fitmaxiter)
|
||||
} else {
|
||||
# use number of kernels determined by baselinefit()
|
||||
eds.fit[[s]] <- pkdecompint(eds.base, intnum = s,
|
||||
k = eds.base$npks[s], maxiter = fitmaxiter)
|
||||
}
|
||||
# Setup the dataframe that makes up the peak table
|
||||
for (kernel in 1:eds.fit[[s]]$num.ker) {
|
||||
eds.fit.parpk <- rbind(eds.fit.parpk,
|
||||
data.frame(peak = factor(eds.fit[[s]]$intnr),
|
||||
kernel = factor(kernel),
|
||||
x = eds.fit[[s]]$parpks[kernel, "loc"],
|
||||
height = eds.fit[[s]]$parpks[kernel, "height"],
|
||||
area = eds.fit[[s]]$parpks[kernel, "intens"],
|
||||
fwhm = eds.fit[[s]]$parpks[kernel, "FWHM"],
|
||||
m = eds.fit[[s]]$parpks[kernel, "m"],
|
||||
accept = eds.fit[[s]]$accept))
|
||||
eds.fit.fitpk <- rbind(eds.fit.fitpk,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = factor(kernel),
|
||||
x = eds.fit[[s]]$x,
|
||||
y = eds.fit[[s]]$fitpk[kernel, ]))
|
||||
}
|
||||
eds.nobasl <- rbind(eds.nobasl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
x = eds.fit[[s]]$x,
|
||||
y = eds.fit[[s]]$y))
|
||||
}
|
||||
|
||||
|
||||
|
||||
return(list(eds.base = eds.base,
|
||||
eds.peaks = eds.pks,
|
||||
eds.fit.parpk = eds.fit.parpk,
|
||||
eds.fit.fitpk = eds.fit.fitpk,
|
||||
eds.nobasl = eds.nobasl))
|
||||
|
||||
}
|
||||
@ -1,6 +0,0 @@
|
||||
# LEO1550.R
|
||||
# Functions to read and manipulate data from the LEO1550 scanning electron microscope
|
||||
# Taha Ahmed, March 2011
|
||||
|
||||
# CONTENTS
|
||||
# >>>> tifftags2df
|
||||
@ -1,121 +0,0 @@
|
||||
##################################################
|
||||
################ tifftags2df #####################
|
||||
##################################################
|
||||
tifftags2df <- function(tiffimage) {
|
||||
## Description:
|
||||
## Extracts all tags from a TIFF image file
|
||||
## using the tiffinfo tool and stores
|
||||
## a selection of the tags in a dataframe.
|
||||
## Usage:
|
||||
## tifftags2df(fulltiffpath)
|
||||
## Arguments:
|
||||
## tiffimage: character string, the full filename
|
||||
## (with path) to one TIFF file.
|
||||
## Value:
|
||||
## A dataframe with three columns:
|
||||
## sampleid, parameter, value.
|
||||
#
|
||||
substrate.id <- strsplit(basename(tiffimage), "\\.")[[1]][1]
|
||||
tifftags <- system(paste("tiffinfo", tiffimage, sep = " "),
|
||||
intern = TRUE, ignore.stderr = TRUE)
|
||||
|
||||
# Clean certain special characters from the tiff tags strings
|
||||
# If these strings are left untreated, they cause all sorts of weird errors later
|
||||
tifftags.clean <-
|
||||
sub("[ ]+$", "", #trim trailing spaces, if any
|
||||
sub("\\r$", "", #remove trailing \r
|
||||
gsub("\\xb9", "", #remove special character
|
||||
gsub("\\xb8", "", #remove special character
|
||||
gsub("\\xb7", "", #remove special character
|
||||
gsub("\\xb6", "", #remove special character
|
||||
gsub("\\xb5", "", #remove special character
|
||||
gsub("\\xb4", "", #remove special character
|
||||
gsub("\\xb3", "", #remove special character
|
||||
gsub("\\xb2", "", #remove special character
|
||||
gsub("\\xb1", "", #remove special character
|
||||
gsub("\\xb0", "", tifftags,
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T),
|
||||
useBytes=T)))
|
||||
|
||||
|
||||
|
||||
# These are the tags we are looking for
|
||||
tags.df <- data.frame(stringsAsFactors = FALSE, substrate.id, matrix(#
|
||||
# Name # REGEXP # splitchar # unit
|
||||
c("EHT", "AP\\_ACTUALKV", " = ", "\\kilo\\volt",
|
||||
"High current", "DP\\_HIGH\\_CURRENT", " = ", "",
|
||||
"WD", "AP\\_WD", " = ", "\\milli\\metre",
|
||||
"Magnification", "AP\\_MAG", " = ", "", # could warrant special treatment
|
||||
"Brightness", "AP\\_BRIGHTNESS", " = ", "\\percent",
|
||||
"Contrast", "AP\\_CONTRAST", " = ", "\\percent",
|
||||
"Signal A", "DP\\_DETECTOR\\_CHANNEL", " = ", "",
|
||||
"SCM status", "DP\\_SCM\\_STATUS", " = ", "",
|
||||
"Specimen current", "AP\\_SCM", " = ", "\\femto\\ampere",
|
||||
# Beam
|
||||
#"Filament curent", "AP\\_ACTUALCURRENT", " = ", "\\ampere",
|
||||
"Stigmation Y", "AP\\_STIG\\_Y", " = ", "\\percent",
|
||||
"Stigmation X", "AP\\_STIG\\_X", " = ", "\\percent",
|
||||
"Aperture align Y", "AP\\_APERTURE\\_ALIGN\\_Y", " = ", "\\percent",
|
||||
"Aperture align X", "AP\\_APERTURE\\_ALIGN\\_X", " = ", "\\percent",
|
||||
"Aperture size", "AP\\_APERTURESIZE", " = ", "\\micro\\metre",
|
||||
"Beam shift Y", "AP\\_BEAMSHIFT\\_Y", " = ", "\\percent",
|
||||
"Beam shift X", "AP\\_BEAMSHIFT\\_X", " = ", "\\percent",
|
||||
#"Beam offset Y", "AP\\_BEAM\\_OFFSET\\_Y", " = ", "\\nano\\metre",
|
||||
#"Beam offset X", "AP\\_BEAM\\_OFFSET\\_X", " = ", "\\nano\\metre",
|
||||
# Stage
|
||||
"Track Z", "DP\\_TRACK\\_Z", " = ", "",
|
||||
"Stage at Z", "AP\\_STAGE\\_AT\\_Z", " = ", "\\milli\\metre",
|
||||
"Stage at Y", "AP\\_STAGE\\_AT\\_Y", " = ", "\\milli\\metre",
|
||||
"Stage at X", "AP\\_STAGE\\_AT\\_X", " = ", "\\milli\\metre",
|
||||
# Tilt
|
||||
"Stage tilted?", "DP\\_STAGE\\_TILTED", " = ", "",
|
||||
"Tilt angle", "AP\\_TILT\\_ANGLE", " = ", "",
|
||||
"Tilt axis", "AP\\_TILT\\_AXIS", " = ", "",
|
||||
# Image
|
||||
"Scan speed", "DP\\_SCANRATE", " = ", "",
|
||||
"Cycle time", "AP\\_FRAME\\_TIME", " = ", "\\second",
|
||||
"Freeze on", "DP\\_FREEZE\\_ON", " = ", "",
|
||||
"Dwell time", "DP\\_DWELL\\_TIME", " = ", "\\nano\\second",
|
||||
"Noise reduction", "DP\\_NOISE\\_REDUCTION", " = ", "",
|
||||
"Frames to integrate", "AP\\_FRAME\\_INT\\_COUNT", " = ", "",
|
||||
"Frames to average", "AP\\_FRAME\\_AVERAGE\\_COUNT", " = ", "",
|
||||
#"Profile width", "AP\\_PROFILE\\_W", " = ", "\\micro\\metre",
|
||||
"Pixel size", "AP\\_PIXEL\\_SIZE", " = ", "\\nano\\metre",
|
||||
# System
|
||||
"Gun vacuum", "AP\\_COLUMN\\_VAC", " = ", "\\milli\\bar",
|
||||
"System vacuum", "AP\\_SYSTEM\\_VAC", " = ", "\\milli\\bar",
|
||||
"Filament age", "AP\\_FILAMENT\\_AGE", " = ", "\\hour",
|
||||
# Misc
|
||||
"Photo no.", "AP\\_PHOTO\\_NUMBER", " = ", "",
|
||||
"Date", "AP\\_DATE", " :", "",
|
||||
"Time", "AP\\_TIME", " :", ""),
|
||||
ncol = 4, byrow = T))
|
||||
names(tags.df) <- c("sampleid", "name", "regexp", "splitchar", "unit")
|
||||
|
||||
|
||||
for (i in 1:dim(tags.df)[1]) {
|
||||
current.tag <- which(regexpr(tags.df$regexp[i], tifftags.clean) == 1) + 1
|
||||
value.tmp <- strsplit(tifftags.clean[current.tag], split = tags.df$splitchar[i])[[1]][2]
|
||||
# Remove leading spaces from tags.df$value
|
||||
tags.df$value[i] <- sub("^\\s+", "", value.tmp, useBytes = T)
|
||||
if (tags.df$unit[i] != "") {
|
||||
tags.df$value[i] <- paste("\\SI{", strsplit(tags.df$value[i], split = " ")[[1]][1], "}{", tags.df$unit[i], "}", sep = "")
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
tags <- data.frame(stringsAsFactors = FALSE,
|
||||
sampleid = tags.df$sampleid,
|
||||
parameter = tags.df$name,
|
||||
value = tags.df$value)
|
||||
|
||||
return(tags)
|
||||
}
|
||||
@ -1,157 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
source("/home/taha/chepec/chetex/common/R/common/int2padstr.R")
|
||||
|
||||
|
||||
##################################################
|
||||
#################### OO2df #######################
|
||||
##################################################
|
||||
OO2df <- function(datafile) {
|
||||
## Description:
|
||||
##
|
||||
##
|
||||
##
|
||||
##
|
||||
## Usage:
|
||||
## OO2df(datafile)
|
||||
## Arguments:
|
||||
## datafile: text string with full path to TXT file
|
||||
## containing single or multiple data ranges
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ wavelength : num
|
||||
## $ counts : num
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $
|
||||
## $ cps ?
|
||||
#
|
||||
# range.header.start.rexp <- "^; \\(Data for Range" #regexp
|
||||
# range.header.end.rexp <- "^_2THETA[^=]" #regexp
|
||||
|
||||
range.data.start.rexp <- ">+Begin[\\s\\w]*<+"
|
||||
range.data.end.rexp <- ">+End[\\s\\w]*<+"
|
||||
|
||||
# Read the input file
|
||||
dfile <- file(datafile, "r")
|
||||
# Note that readLines apparently completely skips empty lines.
|
||||
# That causes line numbers to not match between source and f vector.
|
||||
f <- readLines(dfile, n=-1) # read _all_ lines from data file
|
||||
close(dfile)
|
||||
|
||||
# Fetch a sampleid for the current job
|
||||
sampleid <- ProvideSampleId(datafile)
|
||||
|
||||
# # Look for header start rows
|
||||
# range.header.start.rows <- which(regexpr(range.header.start.rexp, f) == 1)
|
||||
# # Look for header end rows
|
||||
# range.header.end.rows <- which(regexpr(range.header.end.rexp, f) == 1)
|
||||
|
||||
# Look for data start marker line
|
||||
range.data.start.rows <- which(regexpr(range.data.start.rexp, f, perl = TRUE) == 1) + 1
|
||||
# Look for data end marker line
|
||||
range.data.end.rows <- which(regexpr(range.data.end.rexp, f, perl = TRUE) == 1) - 1
|
||||
|
||||
# Calculate number of ranges
|
||||
ranges.total <-
|
||||
ifelse(length(range.data.start.rows) == length(range.data.end.rows),
|
||||
length(range.data.start.rows),
|
||||
NA) #why would they not be equal?
|
||||
if (is.na(ranges.total)) {
|
||||
# Obviously something bad happened.
|
||||
# Do something about it. echo an error message perhaps.
|
||||
# But why would they not be equal?
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Header always precedes start of data
|
||||
range.header.end.rows <- range.data.start.rows - 2
|
||||
if (ranges.total > 1) {
|
||||
range.header.start.rows <- c(1, range.data.end.rows[2:length(range.data.end.rows)])
|
||||
} else {
|
||||
# Data in fact only contains one range
|
||||
range.header.start.rows <- 1
|
||||
}
|
||||
|
||||
# Extract headers (as-is) and put them in a list (by range)
|
||||
headers.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
headers.raw[[range]] <- f[range.header.start.rows[range]:range.header.end.rows[range]]
|
||||
}
|
||||
|
||||
####
|
||||
|
||||
# Extract data (as-is) and put it an list (by range)
|
||||
data.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
data.raw[[range]] <- f[range.data.start.rows[range]:range.data.end.rows[range]]
|
||||
# Replace commas by dots
|
||||
data.raw[[range]] <- gsub(",", "\\.", data.raw[[range]])
|
||||
}
|
||||
|
||||
|
||||
# Specify header parameters to include in dataframe
|
||||
header.param.rexp <-
|
||||
c(DateTime = "^Date:",
|
||||
IntegrationTime = "^Integration Time \\(usec\\):",
|
||||
n_Averaged = "^Spectra Averaged:",
|
||||
Boxcar = "^Boxcar Smoothing:",
|
||||
CorrElectricDark = "^Correct for Electrical Dark:",
|
||||
StrobeLampEnabled = "^Strobe/Lamp Enabled:",
|
||||
CorrDetectorNonLin = "^Correct for Detector Non-linearity:",
|
||||
CorrStrayLight = "^Correct for Stray Light:",
|
||||
n_Pixels = "^Number of Pixels")
|
||||
|
||||
# Collect data and header parameters in dataframes, by range in a list
|
||||
data <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
zz <- textConnection(data.raw[[range]], "r")
|
||||
data[[range]] <- data.frame(stringsAsFactors = FALSE,
|
||||
sampleid,
|
||||
int2padstr(range, "0", 3),
|
||||
matrix(scan(zz, what = character()), ncol = 2, byrow = T))
|
||||
close(zz)
|
||||
# Collect header parameters
|
||||
for (param in 1:length(header.param.rexp)) {
|
||||
data[[range]] <-
|
||||
cbind(stringsAsFactors = FALSE,
|
||||
data[[range]],
|
||||
# Matches any word, digit, plus, or minus character
|
||||
# surrounded by parentheses at the end of the string
|
||||
sub("\\s\\([\\w\\d\\+\\-]+\\)$", "",
|
||||
strsplit(headers.raw[[range]][which(regexpr(unname(header.param.rexp[param]),
|
||||
headers.raw[[range]]) == 1)], ": ")[[1]][2],
|
||||
perl = TRUE))
|
||||
}
|
||||
names(data[[range]]) <-
|
||||
c("sampleid", "range", "wavelength", "intensity", names(header.param.rexp))
|
||||
}
|
||||
|
||||
# Create a unified dataframe
|
||||
data.df <- data[[1]]
|
||||
if (ranges.total > 1) {
|
||||
for (range in 2:ranges.total) {
|
||||
data.df <- rbind(data.df, data[[range]])
|
||||
}
|
||||
}
|
||||
|
||||
# Convert the DateTime column to more legibly (and compact) format
|
||||
data.df$DateTime <-
|
||||
format(as.POSIXct(gsub("\\s[A-Z]{4}\\s", " ", data.df$Date),
|
||||
format = "%a %b %d %H:%M:%S %Y"),
|
||||
format = "%Y-%m-%d %H:%M:%S")
|
||||
# Convert wavelength and intensity to numeric format
|
||||
mode(data.df$wavelength) <- "numeric"
|
||||
mode(data.df$intensity) <- "numeric"
|
||||
|
||||
return(data.df)
|
||||
}
|
||||
@ -1,49 +0,0 @@
|
||||
##################################################
|
||||
#################### SP2df #######################
|
||||
##################################################
|
||||
SP2df <- function(datafile) {
|
||||
## Description:
|
||||
## For now just extracting the bare minimum (the data itself) from SP (ASCII) spectra files.
|
||||
## Usage:
|
||||
## SP2df(datafile)
|
||||
## Arguments:
|
||||
## datafile: text string with full path to TXT file
|
||||
## containing single or multiple data ranges
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ wavelength : num
|
||||
## $ intensity : num
|
||||
#
|
||||
range.data.start.rexp <- "\\#DATA"
|
||||
#range.data.end.rexp <- ">+End[\\s\\w]*<+"
|
||||
|
||||
# Read the input file
|
||||
dfile <- file(datafile, "r")
|
||||
# Note that readLines apparently completely skips empty lines.
|
||||
# That causes line numbers to not match between source and f vector.
|
||||
f <- readLines(dfile, n=-1) # read _all_ lines from data file
|
||||
close(dfile)
|
||||
|
||||
# Create a sampleid for the current job (use the folder name)
|
||||
sampleid <- basename(dirname(datafile))
|
||||
|
||||
# Look for data start marker line
|
||||
range.data.start.row <- grep(range.data.start.rexp, f, perl = TRUE) + 1
|
||||
# Data ends one line before EOF
|
||||
range.data.end.row <- length(f) - 1
|
||||
|
||||
# Extract data (as-is)
|
||||
data.raw <- f[range.data.start.row:range.data.end.row]
|
||||
|
||||
|
||||
# Collect data into dataframe
|
||||
zz <- textConnection(data.raw, "r")
|
||||
data <- data.frame(stringsAsFactors = FALSE,
|
||||
sampleid,
|
||||
matrix(scan(zz, what = numeric()), ncol = 2, byrow = T))
|
||||
close(zz)
|
||||
names(data) <- c("sampleid", "wavelength", "intensity")
|
||||
|
||||
return(data)
|
||||
}
|
||||
@ -1,55 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
Raman2df <- function(datafilename) {
|
||||
## Description:
|
||||
## Reads Raman data in ASCII format
|
||||
## (wavenumber, counts)
|
||||
## and returns a dataframe with the original data,
|
||||
## as well as interpolated wavenumber and counts values
|
||||
## (interpolated data is evenly spaced along x-axis)
|
||||
## Usage:
|
||||
## Raman2df(datafilename)
|
||||
## Arguments:
|
||||
## datafilename: text string with full path to experimental file
|
||||
## Value:
|
||||
## Dataframe with the following columns (and no extra attributes):
|
||||
## $ sampleid : chr (id)
|
||||
## $ wavenum : num (measure)
|
||||
## $ counts : num (measure)
|
||||
## $ wnum.interp : num (measure)
|
||||
## $ cts.interp : num (measure)
|
||||
## Note:
|
||||
##
|
||||
#
|
||||
datafile <- file(datafilename, "r")
|
||||
chifile <- readLines(datafile, n = -1) #read all lines of input file
|
||||
close(datafile)
|
||||
#
|
||||
#####
|
||||
sampleid <- ProvideSampleId(datafilename)
|
||||
#
|
||||
ff <- data.frame(NULL)
|
||||
zz <- textConnection(chifile, "r")
|
||||
ff <- rbind(ff, data.frame(stringsAsFactors = FALSE,
|
||||
sampleid,
|
||||
matrix(scan(zz, what = numeric(), sep = "\t"),
|
||||
ncol = 2, byrow = T)))
|
||||
close(zz)
|
||||
names(ff) <- c("sampleid", "wavenum", "counts")
|
||||
# Sort dataframe by increasing wavenumbers
|
||||
ff <- ff[order(ff$wavenum), ]
|
||||
# ... sort the rownames as well
|
||||
row.names(ff) <- seq(1, dim(ff)[1])
|
||||
# Add interpolated, evenly spaced data to dataframe
|
||||
ff <- cbind(ff,
|
||||
wnum.interp = approx(x = ff$wavenum,
|
||||
y = ff$counts,
|
||||
method = "linear",
|
||||
n = length(ff$wavenum))$x,
|
||||
cts.interp = approx(x = ff$wavenum,
|
||||
y = ff$counts,
|
||||
method = "linear",
|
||||
n = length(ff$wavenum))$y)
|
||||
##
|
||||
return(ff)
|
||||
}
|
||||
@ -1,64 +0,0 @@
|
||||
RamanWrapper <-
|
||||
function(data.exp, run, override = FALSE,
|
||||
kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1) {
|
||||
# the override flag is currently not used
|
||||
|
||||
print("... Started RamanWrapper")
|
||||
|
||||
# check if Ramanpk has already completed successfully for the current job
|
||||
current.dirname <- getwd()
|
||||
print(current.dirname)
|
||||
current.filename <- "raman-peak-data.rda"
|
||||
ramandatafile <- paste(current.dirname, current.filename, sep = "/")
|
||||
|
||||
|
||||
|
||||
if (file.exists(ramandatafile)) {
|
||||
print("... Started if-clause 1")
|
||||
|
||||
# File already exists
|
||||
# return the data using load() or data()
|
||||
|
||||
load(file = ramandatafile)
|
||||
|
||||
if (run > length(ramres)) {
|
||||
|
||||
print("... Started if-clause 1:1")
|
||||
|
||||
# the it does not really exist
|
||||
ramres[[run]] <- Ramanpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor)
|
||||
save(ramres, file = ramandatafile)
|
||||
|
||||
print("... Ended if-clause 1:1")
|
||||
}
|
||||
|
||||
print("... Ended if-clause 1")
|
||||
|
||||
return(ramres)
|
||||
} else {
|
||||
|
||||
print("... Started else-clause 1")
|
||||
|
||||
if (!exists("ramres")) {
|
||||
ramres <- list()
|
||||
print("... ramres list created")
|
||||
}
|
||||
|
||||
# Need to call Ramanpk() and save its results to file as above
|
||||
ramres[[run]] <- Ramanpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor)
|
||||
save(ramres, file = ramandatafile)
|
||||
|
||||
print("... Ended else-clause 1")
|
||||
|
||||
return(ramres)
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,4 +0,0 @@
|
||||
I wanted to run the \texttt{diffractometry} peak-analysis code (which can be both time-consuming and prone to stray errors halting compilation) without the need to re-run the peak analysis every time the document is re-compiled.
|
||||
|
||||
For that purpose, I created this wrapper function.
|
||||
Its job is to call \Rfun{Ramanpk()} only if necessary (i.e., peak analysis has not already been performed). That is made possible by saving the results of a successful analysis to a \texttt{raman-peak-data.rda} file in the directory of the report.
|
||||
@ -1,104 +0,0 @@
|
||||
Ramanpk <-
|
||||
function(data.exp, kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1) {
|
||||
|
||||
print("... Starting baseline fitting")
|
||||
|
||||
data.basl <- baselinefit(data.exp,
|
||||
tau = 2.0,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = 200)
|
||||
|
||||
print("... Ended baseline fitting")
|
||||
|
||||
# This loop deals with the output from baselinefit()
|
||||
# It makes a "melted" dataframe in long form for each
|
||||
# separated peak for some baseline parameters
|
||||
data.pks <- data.frame()
|
||||
data.pks.basl <- data.frame()
|
||||
data.pks.pmg <- data.frame()
|
||||
data.pks.spl <- data.frame()
|
||||
peaks <- 1:length(data.basl$npks)
|
||||
|
||||
|
||||
for (s in peaks) {
|
||||
# recorded data in long form by separated peak
|
||||
data.pks <- rbind(data.pks, # column names assigned after loop
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], ]))
|
||||
# the calculated baseline in long form by separated peak
|
||||
data.pks.basl <- rbind(data.pks.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$baseline$basisl[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
# the taut string estimation in long form by separated peak
|
||||
data.pks.pmg <- rbind(data.pks.pmg,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$pmg$fn[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
# the weighted smoothed spline in long form by separated peak
|
||||
data.pks.spl <- rbind(data.pks.spl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$spl$reg[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Column names assigned to d.pks
|
||||
names(data.pks) <- c("peak", "kernel", "x", "y")
|
||||
|
||||
|
||||
# This loop calls pkdecompint() on each peak separately
|
||||
# It makes a "melted" dataframe in long form for:
|
||||
data.fit <- list() # holds pkdecompint output
|
||||
data.fit.fitpk <- data.frame() # contains fitting curves
|
||||
data.fit.parpk <- data.frame() # physical parameters by peak and kernel
|
||||
data.fit.basl <- data.frame() # data with baseline removed
|
||||
peaks <- 1:length(data.basl$npks)
|
||||
for (s in peaks) {
|
||||
######## PKDECOMPINT ########
|
||||
if (length(kerpk) > 1) {
|
||||
# set number of kernels per peak manually
|
||||
data.fit[[s]] <- pkdecompint(data.basl, intnum = s,
|
||||
k = kerpk[s], maxiter = fitmaxiter)
|
||||
} else {
|
||||
# use number of kernels determined by baselinefit()
|
||||
data.fit[[s]] <- pkdecompint(data.basl, intnum = s,
|
||||
k = data.basl$npks[s], maxiter = fitmaxiter)
|
||||
}
|
||||
# Setup the dataframe that makes up the peak table
|
||||
for (kernel in 1:data.fit[[s]]$num.ker) {
|
||||
data.fit.parpk <- rbind(data.fit.parpk,
|
||||
data.frame(peak = factor(data.fit[[s]]$intnr),
|
||||
kernel = factor(kernel),
|
||||
x = data.fit[[s]]$parpks[kernel, "loc"],
|
||||
height = data.fit[[s]]$parpks[kernel, "height"],
|
||||
area = data.fit[[s]]$parpks[kernel, "intens"],
|
||||
fwhm = data.fit[[s]]$parpks[kernel, "FWHM"],
|
||||
m = data.fit[[s]]$parpks[kernel, "m"],
|
||||
accept = data.fit[[s]]$accept))
|
||||
data.fit.fitpk <- rbind(data.fit.fitpk,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = factor(kernel),
|
||||
x = data.fit[[s]]$x,
|
||||
y = data.fit[[s]]$fitpk[kernel, ]))
|
||||
}
|
||||
data.fit.basl <- rbind(data.fit.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
x = data.fit[[s]]$x,
|
||||
y = data.fit[[s]]$y))
|
||||
}
|
||||
|
||||
|
||||
|
||||
return(list(data.basl = data.basl,
|
||||
data.peaks = data.pks,
|
||||
data.fit.parpk = data.fit.parpk,
|
||||
data.fit.fitpk = data.fit.fitpk,
|
||||
data.fit.basl = data.fit.basl))
|
||||
}
|
||||
@ -1,7 +0,0 @@
|
||||
# Renishaw.R
|
||||
# Functions to read data from the Renishaw Raman spectrometer
|
||||
# Taha Ahmed, Feb 2011
|
||||
|
||||
# CONTENTS
|
||||
# >>>> Raman2df
|
||||
# >>>> Ramanpk
|
||||
@ -1,65 +0,0 @@
|
||||
SiliconWrapper <-
|
||||
function(data.exp, run, override = FALSE,
|
||||
kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1) {
|
||||
# the override flag is currently not used
|
||||
|
||||
print("... Started SiliconWrapper")
|
||||
|
||||
# check if Ramanpk has already completed successfully for the current job
|
||||
current.dirname <- getwd()
|
||||
print(current.dirname)
|
||||
current.filename <- "silicon-peak-data.rda"
|
||||
ramandatafile <- paste(current.dirname, current.filename, sep = "/")
|
||||
|
||||
|
||||
|
||||
if (file.exists(ramandatafile) && !override) {
|
||||
print("... Started if-clause 1")
|
||||
|
||||
# File already exists
|
||||
# return the data using load() or data()
|
||||
|
||||
load(file = ramandatafile)
|
||||
|
||||
if (run > length(ramres)) {
|
||||
|
||||
print("... Started if-clause 1:1")
|
||||
|
||||
# the it does not really exist
|
||||
ramres[[run]] <- Ramanpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor)
|
||||
save(ramres, file = ramandatafile)
|
||||
|
||||
print("... Ended if-clause 1:1")
|
||||
}
|
||||
|
||||
print("... Ended if-clause 1")
|
||||
|
||||
return(ramres)
|
||||
} else {
|
||||
# File does not exist, or override requested
|
||||
|
||||
print("... Started else-clause 1")
|
||||
|
||||
if (!exists("ramres")) {
|
||||
ramres <- list()
|
||||
print("... ramres list created")
|
||||
}
|
||||
|
||||
# Need to call Ramanpk() and save its results to file as above
|
||||
ramres[[run]] <- Ramanpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor)
|
||||
save(ramres, file = ramandatafile)
|
||||
|
||||
print("... Ended else-clause 1")
|
||||
|
||||
return(ramres)
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,21 +0,0 @@
|
||||
# Collection of R functions for analysis of thin-film X-ray diffraction patterns
|
||||
# Maintainer: Taha Ahmed
|
||||
|
||||
# CONTENTS
|
||||
# >>>> matchpdf
|
||||
# >>>> scherrer
|
||||
# >>>> pdf2df
|
||||
# >>>> uxd2mtx
|
||||
# >>>> uxd2df
|
||||
# >>>> muxd2df
|
||||
# >>>> muxd2mtx
|
||||
# >>>> muxd2ls
|
||||
# - REPAIR SHOP
|
||||
# - EliminateKa2
|
||||
# - print.xtable.booktabs
|
||||
# - split.muxd
|
||||
# - strip.ka2
|
||||
# - pearson.beta
|
||||
|
||||
## All XRD reading functions need to be re-written to identify the data instead of identifying metadata
|
||||
## Take inspiration from the functions in CHI.R
|
||||
@ -1,80 +0,0 @@
|
||||
##################################################
|
||||
################## matchpdf ######################
|
||||
##################################################
|
||||
matchpdf <- function(expcol, pdfrow) {
|
||||
# Function for matching each element of two numeric vectors
|
||||
# to the closest-lying element (numerically) in the other
|
||||
# vector. For example for matching 2th values of PDF to
|
||||
# 2th values of recorded diffractogram.
|
||||
## Args: two vectors to be matched (any special considerations on vector lengths?)
|
||||
## Values: A list of 5 numeric vectors. See description inside function.
|
||||
# A word of caution. This function creates a matrix with dimensions pdfrow * expcol.
|
||||
# Usually both of these vectors are rather short. But watch out if large vectors
|
||||
# are submitted --- the looping below could make this function very slow.
|
||||
#
|
||||
diff.thth <- matrix(NA, length(pdfrow), length(expcol))
|
||||
diff.indx <- matrix(NA, length(pdfrow), length(expcol))
|
||||
|
||||
for (pdf.row in 1:length(pdfrow)) {
|
||||
# For every row, calculate differences between that pdf value and all exp values
|
||||
diff.thth[pdf.row, ] <- expcol - pdfrow[pdf.row]
|
||||
}
|
||||
#
|
||||
# Now we go on to find the minimum along each row of diff.thth
|
||||
# and set all other values to some arbitrary value
|
||||
diff.rmin <- diff.thth
|
||||
for (pdf.row in 1:dim(diff.thth)[1]) {
|
||||
for (exp.col in 1:dim(diff.thth)[2]) {
|
||||
if (abs(diff.thth[pdf.row, exp.col]) != min(abs(diff.thth[pdf.row, ]))) {
|
||||
diff.rmin[pdf.row, exp.col] <- Inf
|
||||
}
|
||||
}
|
||||
}
|
||||
#
|
||||
# We likewise find the minimum along each column and set any other,
|
||||
# non-Inf values to Inf
|
||||
diff.cmin <- diff.rmin
|
||||
for (exp.col in 1:dim(diff.rmin)[2]) {
|
||||
for (pdf.row in 1:dim(diff.rmin)[1]) {
|
||||
if (abs(diff.rmin[pdf.row, exp.col]) != min(abs(diff.rmin[, exp.col]))) {
|
||||
diff.cmin[pdf.row, exp.col] <- Inf
|
||||
}
|
||||
}
|
||||
}
|
||||
#
|
||||
# The matrix now contains at most one non-Inf value per row and column.
|
||||
# To make the use of colSums and rowSums possible (next step) all Inf are set zero,
|
||||
# and all matches are set to 1.
|
||||
for (pdf.row in 1:dim(diff.cmin)[1]) {
|
||||
for (exp.col in 1:dim(diff.cmin)[2]) {
|
||||
if ((diff.cmin[pdf.row, exp.col]) != Inf) {
|
||||
diff.indx[pdf.row, exp.col] <- 1
|
||||
} else {
|
||||
diff.indx[pdf.row, exp.col] <- 0
|
||||
}
|
||||
}
|
||||
}
|
||||
#
|
||||
# Attach error message to return list if sum(rowSums(diff.indx)) == sum(colSums(diff.indx))
|
||||
matchpdf.err.msg <-"matchpdf() failed to complete. It seems rowSums != colSums."
|
||||
mtch <- matchpdf.err.msg
|
||||
#
|
||||
# Check that sum(rowSums(diff.indx)) == sum(colSums(diff.indx))
|
||||
if (sum(rowSums(diff.indx)) == sum(colSums(diff.indx))) {
|
||||
# Reset mtch
|
||||
mtch <- list()
|
||||
mtch <- list(csums = colSums(diff.indx),
|
||||
rsums = rowSums(diff.indx),
|
||||
expthth = expcol[colSums(diff.indx) != 0],
|
||||
pdfthth = pdfrow[rowSums(diff.indx) != 0],
|
||||
deltathth = expcol[colSums(diff.indx) != 0] - pdfrow[rowSums(diff.indx) != 0])
|
||||
# List of 5
|
||||
# $ csums : num - consisting of ones and zeroes. Shows you which positions of expcol matched.
|
||||
# $ rsums : num - consisting of ones and zeroes. Shows you which positions of pdfrow matched.
|
||||
# $ expthth : num - consisting of the matching elements of expcol.
|
||||
# $ pdfthth : num - consisting of the matching elements of pdfrow.
|
||||
# $ deltathth : num - element-wise difference of expcol and pdfrow.
|
||||
}
|
||||
#
|
||||
return(mtch)
|
||||
}
|
||||
@ -1,154 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
source("/home/taha/chepec/chetex/common/R/common/int2padstr.R")
|
||||
|
||||
|
||||
##################################################
|
||||
#!!!!!!!!!!!!!!!!!# mraw2df #!!!!!!!!!!!!!!!!!!!!#
|
||||
##################################################
|
||||
mraw2df <- function(txtfile) {
|
||||
## Description:
|
||||
## Reads TXT files with one or multiple ranges (XRD commander "save as txt")
|
||||
## Extracts both data (thth, intensity) and parameters
|
||||
## Usage:
|
||||
## mraw2df(txtfile)
|
||||
## Arguments:
|
||||
## txtfile: text string with full path to TXT file, which may
|
||||
## containing single or multiple data ranges
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ thth : num
|
||||
## $ counts (or cps) : num
|
||||
## $ steptime : num
|
||||
## $ stepsize : num
|
||||
## $ theta : num
|
||||
## $ khi : num
|
||||
## $ phi : num
|
||||
## $ x : num
|
||||
## $ y : num
|
||||
## $ z : num
|
||||
## $ divergence : num
|
||||
## $ antiscatter : num
|
||||
## $ cps (or counts) : num
|
||||
#
|
||||
range.header.start.rexp <- "^; \\(Data for Range" #regexp
|
||||
range.header.end.rexp <- "^_2THETA[^=]" #regexp
|
||||
|
||||
# Read the input multirange file
|
||||
ufile <- file(uxdfile, "r")
|
||||
# Note that readLines apparently completely skips empty lines.
|
||||
# In that case line numbers do not match between source and f.
|
||||
f <- readLines(ufile, n=-1) #read _all_ lines from UXD file
|
||||
close(ufile)
|
||||
|
||||
# Fetch a sampleid for the current job
|
||||
sampleid <- ProvideSampleId(uxdfile)
|
||||
|
||||
# Look for header start rows
|
||||
range.header.start.rows <- which(regexpr(range.header.start.rexp, f) == 1)
|
||||
# Look for header end rows
|
||||
range.header.end.rows <- which(regexpr(range.header.end.rexp, f) == 1)
|
||||
|
||||
# Calculate number of ranges
|
||||
ranges.total <-
|
||||
ifelse(length(range.header.start.rows) == length(range.header.end.rows),
|
||||
length(range.header.start.rows),
|
||||
NA) #why would they not be equal?
|
||||
if (is.na(ranges.total)) {
|
||||
# Obviously something bad happened.
|
||||
# Do something about it. echo an error message perhaps.
|
||||
# But why would they not be equal?
|
||||
|
||||
}
|
||||
|
||||
# Determine whether we have COUNTS or COUNTS PER SECOND in current UXD-file
|
||||
# Assuming it is the same for all ranges in this job (a safe assumption).
|
||||
if (f[range.header.end.rows][1] == "_2THETACOUNTS") {
|
||||
# we got counts
|
||||
counts.flag <- TRUE
|
||||
cps.flag <- FALSE
|
||||
}
|
||||
if (f[range.header.end.rows][1] == "_2THETACPS") {
|
||||
# we got counts per second
|
||||
counts.flag <-FALSE
|
||||
cps.flag <- TRUE
|
||||
}
|
||||
|
||||
# Extract headers (as-is) and put them in a list (by range)
|
||||
headers.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
headers.raw[[range]] <- f[range.header.start.rows[range]:range.header.end.rows[range]]
|
||||
}
|
||||
|
||||
# Data always start on the row after header end
|
||||
range.data.start.rows <- range.header.end.rows + 1
|
||||
# Data end rows precedes header with one row, except for the first range
|
||||
# But only if data contained more than one range, obviously. Let's make the code check for that
|
||||
if (ranges.total > 1) {
|
||||
range.data.end.rows <- c(range.header.start.rows[2:length(range.header.start.rows)] - 1, length(f))
|
||||
} else {
|
||||
# Data in fact only contains one range
|
||||
range.data.end.rows <- length(f)
|
||||
}
|
||||
|
||||
####
|
||||
|
||||
# Extract data (as-is) and put it an list (by range)
|
||||
data.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
data.raw[[range]] <- f[range.data.start.rows[range]:range.data.end.rows[range]]
|
||||
}
|
||||
|
||||
# Specify header parameters to include in dataframe
|
||||
header.param.rexp <- c(steptime = "^_STEPTIME=",
|
||||
stepsize = "^_STEPSIZE=",
|
||||
theta = "^_THETA=",
|
||||
khi = "^_KHI=",
|
||||
phi = "^_PHI=",
|
||||
x = "^_X=",
|
||||
y = "^_Y=",
|
||||
z = "^_Z=",
|
||||
divergence = "^_DIVERGENCE=",
|
||||
antiscatter = "^_ANTISCATTER=")
|
||||
|
||||
# Collect data and header parameters in dataframes, by range in a list
|
||||
data <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
zz <- textConnection(data.raw[[range]], "r")
|
||||
data[[range]] <- data.frame(stringsAsFactors = F,
|
||||
sampleid,
|
||||
int2padstr(range, "0", 3),
|
||||
matrix(scan(zz, what = numeric()), ncol = 2, byrow = T))
|
||||
close(zz)
|
||||
# Collect header parameters
|
||||
for (param in 1:length(header.param.rexp)) {
|
||||
data[[range]] <- cbind(data[[range]],
|
||||
as.numeric(strsplit(headers.raw[[range]][which(regexpr(unname(header.param.rexp[param]),
|
||||
headers.raw[[range]]) == 1)], "=")[[1]][2]))
|
||||
}
|
||||
names(data[[range]]) <-
|
||||
c("sampleid", "range", "thth", ifelse(counts.flag, "counts", "cps"), names(header.param.rexp))
|
||||
}
|
||||
|
||||
# Calculate the other of the pair counts <-> cps
|
||||
if (counts.flag) {
|
||||
for (range in 1:ranges.total) {
|
||||
data[[range]] <- cbind(data[[range]], cps = data[[range]]$counts / data[[range]]$steptime)
|
||||
}
|
||||
}
|
||||
if (cps.flag) {
|
||||
for (range in 1:ranges.total) {
|
||||
data[[range]] <- cbind(data[[range]], counts = data[[range]]$cps * data[[range]]$steptime)
|
||||
}
|
||||
}
|
||||
|
||||
# Return a unified dataframe
|
||||
data.df <- data[[1]]
|
||||
if (ranges.total > 1) {
|
||||
for (range in 2:ranges.total) {
|
||||
data.df <- rbind(data.df, data[[range]])
|
||||
}
|
||||
}
|
||||
|
||||
return(data.df)
|
||||
}
|
||||
@ -1,156 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
source("/home/taha/chepec/chetex/common/R/common/int2padstr.R")
|
||||
|
||||
|
||||
##################################################
|
||||
################### muxd2df ######################
|
||||
##################################################
|
||||
muxd2df <- function(uxdfile) {
|
||||
## Description:
|
||||
## Reads UXD files with one or multiple ranges (converted using XCH v1.0)
|
||||
## Extracts both data (thth, intensity) and parameters
|
||||
## Also automatically calculates cps if counts are present, and vice versa
|
||||
## (note that this depends on specific strings in the UXD format).
|
||||
## Usage:
|
||||
## muxd2df(uxdfile)
|
||||
## Arguments:
|
||||
## uxdfile: text string with full path to UXD file, which may
|
||||
## containing single or multiple data ranges
|
||||
## Value:
|
||||
## Dataframe with the following columns:
|
||||
## $ sampleid : chr
|
||||
## $ thth : num
|
||||
## $ counts (or cps) : num
|
||||
## $ steptime : num
|
||||
## $ stepsize : num
|
||||
## $ theta : num
|
||||
## $ khi : num
|
||||
## $ phi : num
|
||||
## $ x : num
|
||||
## $ y : num
|
||||
## $ z : num
|
||||
## $ divergence : num
|
||||
## $ antiscatter : num
|
||||
## $ cps (or counts) : num
|
||||
#
|
||||
range.header.start.rexp <- "^; \\(Data for Range" #regexp
|
||||
range.header.end.rexp <- "^_2THETA[^=]" #regexp
|
||||
|
||||
# Read the input multirange file
|
||||
ufile <- file(uxdfile, "r")
|
||||
# Note that readLines apparently completely skips empty lines.
|
||||
# In that case line numbers do not match between source and f.
|
||||
f <- readLines(ufile, n=-1) #read _all_ lines from UXD file
|
||||
close(ufile)
|
||||
|
||||
# Fetch a sampleid for the current job
|
||||
sampleid <- ProvideSampleId(uxdfile)
|
||||
|
||||
# Look for header start rows
|
||||
range.header.start.rows <- which(regexpr(range.header.start.rexp, f) == 1)
|
||||
# Look for header end rows
|
||||
range.header.end.rows <- which(regexpr(range.header.end.rexp, f) == 1)
|
||||
|
||||
# Calculate number of ranges
|
||||
ranges.total <-
|
||||
ifelse(length(range.header.start.rows) == length(range.header.end.rows),
|
||||
length(range.header.start.rows),
|
||||
NA) #why would they not be equal?
|
||||
if (is.na(ranges.total)) {
|
||||
# Obviously something bad happened.
|
||||
# Do something about it. echo an error message perhaps.
|
||||
# But why would they not be equal?
|
||||
|
||||
}
|
||||
|
||||
# Determine whether we have COUNTS or COUNTS PER SECOND in current UXD-file
|
||||
# Assuming it is the same for all ranges in this job (a safe assumption).
|
||||
if (f[range.header.end.rows][1] == "_2THETACOUNTS") {
|
||||
# we got counts
|
||||
counts.flag <- TRUE
|
||||
cps.flag <- FALSE
|
||||
}
|
||||
if (f[range.header.end.rows][1] == "_2THETACPS") {
|
||||
# we got counts per second
|
||||
counts.flag <-FALSE
|
||||
cps.flag <- TRUE
|
||||
}
|
||||
|
||||
# Extract headers (as-is) and put them in a list (by range)
|
||||
headers.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
headers.raw[[range]] <- f[range.header.start.rows[range]:range.header.end.rows[range]]
|
||||
}
|
||||
|
||||
# Data always start on the row after header end
|
||||
range.data.start.rows <- range.header.end.rows + 1
|
||||
# Data end rows precedes header with one row, except for the first range
|
||||
# But only if data contained more than one range, obviously. Let's make the code check for that
|
||||
if (ranges.total > 1) {
|
||||
range.data.end.rows <- c(range.header.start.rows[2:length(range.header.start.rows)] - 1, length(f))
|
||||
} else {
|
||||
# Data in fact only contains one range
|
||||
range.data.end.rows <- length(f)
|
||||
}
|
||||
|
||||
####
|
||||
|
||||
# Extract data (as-is) and put it an list (by range)
|
||||
data.raw <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
data.raw[[range]] <- f[range.data.start.rows[range]:range.data.end.rows[range]]
|
||||
}
|
||||
|
||||
# Specify header parameters to include in dataframe
|
||||
header.param.rexp <- c(steptime = "^_STEPTIME=",
|
||||
stepsize = "^_STEPSIZE=",
|
||||
theta = "^_THETA=",
|
||||
khi = "^_KHI=",
|
||||
phi = "^_PHI=",
|
||||
x = "^_X=",
|
||||
y = "^_Y=",
|
||||
z = "^_Z=",
|
||||
divergence = "^_DIVERGENCE=",
|
||||
antiscatter = "^_ANTISCATTER=")
|
||||
|
||||
# Collect data and header parameters in dataframes, by range in a list
|
||||
data <- list()
|
||||
for (range in 1:ranges.total) {
|
||||
zz <- textConnection(data.raw[[range]], "r")
|
||||
data[[range]] <- data.frame(stringsAsFactors = F,
|
||||
sampleid,
|
||||
int2padstr(range, "0", 3),
|
||||
matrix(scan(zz, what = numeric()), ncol = 2, byrow = T))
|
||||
close(zz)
|
||||
# Collect header parameters
|
||||
for (param in 1:length(header.param.rexp)) {
|
||||
data[[range]] <- cbind(data[[range]],
|
||||
as.numeric(strsplit(headers.raw[[range]][which(regexpr(unname(header.param.rexp[param]),
|
||||
headers.raw[[range]]) == 1)], "=")[[1]][2]))
|
||||
}
|
||||
names(data[[range]]) <-
|
||||
c("sampleid", "range", "thth", ifelse(counts.flag, "counts", "cps"), names(header.param.rexp))
|
||||
}
|
||||
|
||||
# Calculate the other of the pair counts <-> cps
|
||||
if (counts.flag) {
|
||||
for (range in 1:ranges.total) {
|
||||
data[[range]] <- cbind(data[[range]], cps = data[[range]]$counts / data[[range]]$steptime)
|
||||
}
|
||||
}
|
||||
if (cps.flag) {
|
||||
for (range in 1:ranges.total) {
|
||||
data[[range]] <- cbind(data[[range]], counts = data[[range]]$cps * data[[range]]$steptime)
|
||||
}
|
||||
}
|
||||
|
||||
# Return a unified dataframe
|
||||
data.df <- data[[1]]
|
||||
if (ranges.total > 1) {
|
||||
for (range in 2:ranges.total) {
|
||||
data.df <- rbind(data.df, data[[range]])
|
||||
}
|
||||
}
|
||||
|
||||
return(data.df)
|
||||
}
|
||||
@ -1,38 +0,0 @@
|
||||
##################################################
|
||||
################### muxd2ls ######################
|
||||
##################################################
|
||||
muxd2ls <- function(uxdfile) {
|
||||
# Function that reads an UXD file which contains several ranges
|
||||
# (created in a programmed run, for example)
|
||||
# Arguments
|
||||
# :: uxdfile (filename with extension)
|
||||
# Returns: List of matrices, as many as there were ranges
|
||||
# Requires: ??
|
||||
# See extensive comments in muxd2mtx()
|
||||
|
||||
cchar <- "[;_]" #comment characters used in Bruker's UXD
|
||||
cdata <- "[0-9]"
|
||||
|
||||
ufile <- file(uxdfile, "r")
|
||||
f <- readLines(ufile, n=-1) #read _all_ lines from UXD file
|
||||
close(ufile)
|
||||
|
||||
wh <- regexpr(cchar, f)
|
||||
mh <- wh[1:length(wh)]
|
||||
|
||||
i <- seq(1, length(mh) - 1, 1)
|
||||
j <- seq(2, length(mh), 1)
|
||||
starts <- which(mh[i] == 1 & mh[j] != 1) + 1
|
||||
ends <- which(mh[i] != 1 & mh[j] == 1)
|
||||
ends <- c(ends, length(mh))
|
||||
|
||||
ff <- list()
|
||||
for (s in 1:length(starts)) {
|
||||
zz <- textConnection(f[starts[s]:ends[s]], "r")
|
||||
ms <- matrix(scan(zz, what = numeric()), ncol = 2, byrow = T)
|
||||
close(zz)
|
||||
ff[[s]] <- ms
|
||||
}
|
||||
# Return list of matrices
|
||||
return(ff)
|
||||
}
|
||||
@ -1,57 +0,0 @@
|
||||
##################################################
|
||||
################### muxd2mtx #####################
|
||||
##################################################
|
||||
muxd2mtx <- function(uxdfile, range) {
|
||||
# Function that reads an UXD file which contains several ranges
|
||||
# (created in a programmed run, for example)
|
||||
# Arguments
|
||||
# :: uxdfile (filename with extension)
|
||||
# :: range (an integer, the number of ranges in the uxdfile)
|
||||
# Returns: matrix with 2 columns
|
||||
|
||||
cchar <- "[;_]" #regexpr matching the comment characters used in Bruker's UXD
|
||||
cdata <- "[0-9]" #regexpr matching one character of any digit
|
||||
# Create filenames for the output # no longer used, return dataframe instead
|
||||
#datafile <- paste(uxdfile,"-",range.descriptor,".data",sep="")
|
||||
|
||||
# Read the input multirange file
|
||||
ufile <- file(uxdfile, "r")
|
||||
f <- readLines(ufile, n=-1) #read _all_ lines from UXD file
|
||||
close(ufile)
|
||||
|
||||
# We identify header rows. We will assume that all other rows are data
|
||||
wh <- regexpr(cchar, f)
|
||||
# length(wh) equals length(f)
|
||||
# wh has either 1 or -1 for each element
|
||||
# value = 1 means the first character of that row is cchar (row is header)
|
||||
# value =-1 means absence of cchar on that row (row is data)
|
||||
|
||||
# Since wh contains some attributes (given by regexpr function), we strip out everything
|
||||
# but the index vector and assign it to the new vector mh
|
||||
mh <- wh[1:length(wh)] # this gives you the corresponding index vector
|
||||
|
||||
#length(mh[mh == -1]) #total number of datarows in uxdfile
|
||||
#mh[mh > 1 | mh < 0] <- 0 #set all header-rows to zero (just to make things easier)
|
||||
|
||||
# Set counters i and j used in assignment below
|
||||
i <- seq(1, length(mh) - 1, 1)
|
||||
j <- seq(2, length(mh), 1)
|
||||
|
||||
starts <- which(mh[i] == 1 & mh[j] != 1) + 1 #start indices
|
||||
ends <- which(mh[i] != 1 & mh[j] == 1) #end indices, except the last
|
||||
ends <- c(ends, length(mh)) #fixes the last index of ends
|
||||
# note that length of starts (or ends) gives the number of ranges
|
||||
|
||||
ff <- matrix(NA, length(f), 2)
|
||||
for (s in 1:range) {
|
||||
zz <- textConnection(f[starts[s]:ends[s]], "r")
|
||||
ff <- rbind(ff, matrix(scan(zz, what = numeric()), ncol=2, byrow=T))
|
||||
close(zz)
|
||||
}
|
||||
|
||||
# Clean up matrix: remove extra rows
|
||||
ff[apply(ff,1,function(x)any(!is.na(x))), ]
|
||||
|
||||
# Return matrix
|
||||
ff
|
||||
}
|
||||
@ -1,70 +0,0 @@
|
||||
pdf2df <- function(pdffile) {
|
||||
# Function for extracting information from ICDD PDF XML-files
|
||||
# For example the PDF files produced by the PDF database at Angstrom's X-ray lab
|
||||
# NOTE: sometimes intensity values are specified as less than some value.
|
||||
# In those cases, the lt sign will be preserved in the column int.Tex.
|
||||
# The intensity column, on the other hand, is numeric and so strips off the lt sign.
|
||||
# ARGS: pdffile (complete path and filename to PDF file)
|
||||
# VALUE: dataframe with 9 columns:
|
||||
# thth angles (numeric),
|
||||
# d (numeric),
|
||||
# h index (numeric),
|
||||
# k index (numeric),
|
||||
# l index (numeric),
|
||||
# hkl indices (string),
|
||||
# hkl.TeX indices formatted for LaTeX (string),
|
||||
# intensity (numeric),
|
||||
# int.TeX intensity formatted for LaTeX (string),
|
||||
# pdfNumber (string)
|
||||
# attr: This function sets the following attributes:
|
||||
# ApplicationName,
|
||||
# ApplicationVersion,
|
||||
# chemicalformula,
|
||||
# empiricalformula,
|
||||
# wavelength
|
||||
#
|
||||
require(XML)
|
||||
doc <- xmlTreeParse(pdffile)
|
||||
pdf <- xmlRoot(doc)
|
||||
rmchar <- "[^0-9]"
|
||||
#
|
||||
angles <- data.frame(NULL)
|
||||
for (i in 1:length(pdf[["graphs"]][["stick_series"]])) {
|
||||
angles <- rbind(angles, data.frame(stringsAsFactors = FALSE,#
|
||||
thth = as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["theta"]])),
|
||||
d = as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["da"]])),
|
||||
h = as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["h"]])),
|
||||
k = as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["k"]])),
|
||||
l = as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["l"]])),
|
||||
hkl = paste(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["h"]]),
|
||||
xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["k"]]),
|
||||
xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["l"]]), sep = ""),
|
||||
hkl.TeX = paste("\\mbox{$", ifelse(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["h"]])) < 0,
|
||||
paste("\\bar{", abs(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["h"]]))),
|
||||
"}", sep = ""),
|
||||
xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["h"]])),
|
||||
"\\,", ifelse(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["k"]])) < 0,
|
||||
paste("\\bar{", abs(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["k"]]))),
|
||||
"}", sep = ""),
|
||||
xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["k"]])),
|
||||
"\\,", ifelse(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["l"]])) < 0,
|
||||
paste("\\bar{", abs(as.numeric(xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["l"]]))),
|
||||
"}", sep = ""),
|
||||
xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["l"]])),
|
||||
"$}", sep = "", collapse = ""),
|
||||
intensity = as.numeric(gsub(rmchar, "", xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["intensity"]]))),
|
||||
int.TeX = paste("{", xmlValue(pdf[["graphs"]][["stick_series"]][[i]][["intensity"]]), "}", sep = ""),
|
||||
pdfNumber = xmlValue(pdf[["pdf_data"]][["pdf_number"]]),
|
||||
formula = gsub("[ ]", "", xmlValue(pdf[["pdf_data"]][["empirical_formula"]]))
|
||||
))
|
||||
}
|
||||
#
|
||||
attr(angles, "ApplicationName") <- xmlAttrs(pdf)[[1]]
|
||||
attr(angles, "ApplicationVersion") <- xmlAttrs(pdf)[[2]]
|
||||
#attr(angles, "pdfNumber") <- xmlValue(pdf[["pdf_data"]][["pdf_number"]])
|
||||
attr(angles, "chemicalformula") <- gsub("[ ]", "", xmlValue(pdf[["pdf_data"]][["chemical_formula"]]))
|
||||
attr(angles, "empiricalformula") <- gsub("[ ]", "", xmlValue(pdf[["pdf_data"]][["empirical_formula"]]))
|
||||
attr(angles, "wavelength") <- as.numeric(xmlValue(pdf[["graphs"]][["wave_length"]]))
|
||||
# Caution: Do not subset. Subsetting causes all attributes to be lost.
|
||||
return(angles)
|
||||
}
|
||||
@ -1,17 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/as.radians.R")
|
||||
|
||||
##################################################
|
||||
################## scherrer ######################
|
||||
##################################################
|
||||
scherrer <- function(integralbreadth, thth, wavelength = 1.54056, shapeconstant = ((4/3)*(pi/6))^(1/3)) {
|
||||
# Function for calculating crystallite grain size from reflection data
|
||||
# ARGS: integralbreadth - vector with integral breadth of reflections (in degrees)
|
||||
# thth - vector with 2theta values of reflections (in degrees)
|
||||
# wavelength - X-ray wavelength used (default 1.54056 A, Cu Ka)
|
||||
# shapeconstant - Scherrer constant (default spherical, ~0.9)
|
||||
# VALUE: vector with size parameters
|
||||
## REQUIRES: as.radians()
|
||||
D <- (shapeconstant * wavelength) / (as.radians(integralbreadth) * cos(as.radians(thth)))
|
||||
# cos() - angles must be in radians, not degrees!
|
||||
return(D) #units of angstrom
|
||||
}
|
||||
@ -1,215 +0,0 @@
|
||||
##################################################
|
||||
################ EliminateKa2 ####################
|
||||
##################################################
|
||||
EliminateKa2 <- function(thth, intensity, lever = 2) {
|
||||
# Args: 2theta vector, numeric
|
||||
# intensity vector, numeric
|
||||
# parameter pairs, numeric between 1 and 6, default 2
|
||||
# "High-quality a2-free patterns can be obtained in most cases
|
||||
# using five (5) or seven (7) pairs of parameters."
|
||||
# "When the step-width is less than 0.01 degrees and the
|
||||
# 2theta angle is high, a large number of parameter pairs
|
||||
# should be used to get accurate results." {Dong1999a}
|
||||
### 1 2 3 4 5 6 - lever
|
||||
### 3 5 7 9 15 25 - corresponding parameter pairs
|
||||
#
|
||||
##### THIS FUNCTION USES THESE OTHER FUNCTIONS #####
|
||||
# REQUIRES: common.R :: as.radians() - converts degrees to radians
|
||||
# REQUIRES: stats::approx - linear interpolation
|
||||
####################################################
|
||||
|
||||
# For test-purposes, use the following data
|
||||
# /home/taha/chepec/laboratory/XRD/0103-instrumentbroadening/100917Th2ThLong-counts.UXD
|
||||
|
||||
##### STILL UNDER CONSTRUCTION ####
|
||||
##### STILL UNDER CONSTRUCTION ####
|
||||
|
||||
startdatapoint <- 4
|
||||
|
||||
|
||||
# Ka1a Ka1b Ka2a Ka2b
|
||||
CuKa.data <- c(1.534753, 1.540596, 1.541058, 1.544410, 1.544721,
|
||||
3.69, 0.44, 0.60, 0.52, 0.62,
|
||||
1.60, 57.07, 7.64, 25.38, 8.31)
|
||||
CuKa <- data.frame(matrix(CuKa.data, ncol=3, byrow=F))
|
||||
names(CuKa) <- c("lambda", "w", "E")
|
||||
row.names(CuKa) <- c("Satellites", "Ka1a", "Ka1b", "Ka2a", "Ka2b")
|
||||
|
||||
|
||||
# The following lever arm weights are from Dong1999a
|
||||
weights <- list()
|
||||
# Three-bar weights
|
||||
weights[[1]] <- c(0.005134296, 0.491686047, 0.003179657)
|
||||
# Five-bar weights
|
||||
weights[[2]] <- c(0.002614410, 0.011928014, 0.480406967,
|
||||
0.002121807, 0.002928802)
|
||||
# Seven-bar weights
|
||||
weights[[3]] <- c(0.001580069, 0.003463773, 0.015533472, 0.422601977,
|
||||
0.053632977, 0.001572467, 0.001615265)
|
||||
# Nine-bar weights
|
||||
weights[[4]] <- c(0.001138001, 0.00195272, 0.004324464,
|
||||
0.019246541, 0.394175823, 0.079159001,
|
||||
-0.003591547, 0.002505604, 0.001089392)
|
||||
# 15-bar weights
|
||||
weights[[5]] <- c(0.000614225, 0.000810836, 0.001134775, 0.001723265, 0.002968405,
|
||||
0.006433676, 0.02575384, 0.345872599, 0.100578092, 0.014493969,
|
||||
-0.004176171, 0.000678688, 0.001610333, 0.000918077, 0.000585391)
|
||||
# 25-bar weights
|
||||
weights[[6]] <- c(0.000349669, 0.000408044, 0.000484578, 0.000587457, 0.000730087,
|
||||
0.000935685, 0.001247401, 0.001753233, 0.002657209, 0.004531817,
|
||||
0.009591103, 0.034998436, 0.2876498, 0.074964321, 0.065000871,
|
||||
0.016762729, -0.00306221, -0.002717412, -0.000902322, 0.000915701,
|
||||
0.001036484, 0.000808199, 0.000539899, 0.000398896, 0.000330325)
|
||||
|
||||
# The following lever arm lengths are from Dong1999a
|
||||
lengths <- list()
|
||||
# Three-bar lengths
|
||||
lengths[[1]] <- c(0.998506815, 0.997503913, 0.996699460)
|
||||
# Five-bar lengths
|
||||
lengths[[2]] <- c(0.998471576, 0.997935524, 0.997503530,
|
||||
0.997163494, 0.996606519)
|
||||
# Seven-bar lengths
|
||||
lengths[[3]] <- c(0.998563433, 0.998204025, 0.997825027, 0.997522195,
|
||||
0.997297615, 0.996844235, 0.996516288)
|
||||
# Nine-bar lengths
|
||||
lengths[[4]] <- c(0.998609749, 0.998334027, 0.998054914,
|
||||
0.99776062, 0.997527844, 0.997327154,
|
||||
0.997028978, 0.996734639, 0.99646335)
|
||||
# 15-bar lengths
|
||||
lengths[[5]] <- c(0.998671599, 0.99850911, 0.998346447, 0.998183442, 0.998019704,
|
||||
0.997854063, 0.997680649, 0.997533314, 0.997377391, 0.997266106,
|
||||
0.997060614, 0.996888005, 0.996741151, 0.996583672, 0.996418168)
|
||||
# 25-bar lengths
|
||||
lengths[[6]] <- c(0.998706192, 0.998608958, 0.998511721, 0.998414475, 0.998317209,
|
||||
0.998219906, 0.998122538, 0.998025057, 0.997927367, 0.997829244,
|
||||
0.997730044, 0.997626987, 0.997535705, 0.997458223, 0.997346989,
|
||||
0.997277763, 0.997161452, 0.997057942, 0.996982688, 0.99686108,
|
||||
0.996769728, 0.996675255, 0.996578407, 0.996480641, 0.99638324)
|
||||
|
||||
### --- Arguments check
|
||||
# Check that "lever" argument is within bounds
|
||||
if (lever > length(weights) || lever < 1) {
|
||||
# if not, fall back to the default value
|
||||
lever <- 2 # corresponds to 5 parameter pairs
|
||||
}
|
||||
|
||||
### --- Arguments check
|
||||
# Check that vectors are of the same length
|
||||
if (!(length(thth) == length(intensity))) {
|
||||
# If not the same length, abort with error message
|
||||
stop("Arguments thth and intensity have different lengths!")
|
||||
}
|
||||
|
||||
### --- Arguments check
|
||||
if (any(thth <= 0)) {
|
||||
stop("thth vector contains values less-than or equal to zero")
|
||||
}
|
||||
|
||||
### --- Arguments check
|
||||
if (any(intensity < 0)) {
|
||||
stop("intensity vector contains values less than zero")
|
||||
}
|
||||
|
||||
# THIS IS NECESSARY, but overlooked for the moment
|
||||
#int.p.start <- (1 / (1 + (CuKa["Ka2a", "E"] + CuKa["Ka2b", "E"]) /
|
||||
#(CuKa["Ka1a", "E"] + CuKa["Ka1b", "E"]))) * intensity[1:startdatapoint-1]
|
||||
|
||||
# Redefine everything
|
||||
#thth <- thth[startdatapoint:length(thth)]
|
||||
#intensity <- intensity[startdatapoint:length(intensity)]
|
||||
|
||||
|
||||
# Convert from 2theta to theta scale for use in first step of calculations
|
||||
theta <- thth / 2
|
||||
sintheta <- sin(as.radians(theta))
|
||||
|
||||
# Corresponds to equation 10 in Dong1999a
|
||||
# This is based on the assumption that we are supposed to get delta-thth values
|
||||
delta.thth.a <- matrix(0, length(theta), length(weights[[lever]]))
|
||||
for (j in 1:length(lengths[[lever]])) {
|
||||
delta.thth.a[,j] <- 2 * asin(as.radians((lengths[[lever]][j] * sintheta)))
|
||||
}
|
||||
|
||||
# Add the calculated deltas to the recorded thth values
|
||||
# Corresponds to equation 10 in Dong1999a
|
||||
thth.a <- matrix(NA, dim(delta.thth.a)[1], dim(delta.thth.a)[2] + 1)
|
||||
# Flip the delta.thth.a matrix vertically (just for convenience)
|
||||
delta.thth.a <- delta.thth.a[,dim(delta.thth.a)[2]:1]
|
||||
for (j in 2:dim(thth.a)[2]) {
|
||||
thth.a[,1] <- thth
|
||||
thth.a[,j] <- thth + delta.thth.a[,j - 1]
|
||||
}
|
||||
|
||||
# Intensities with interpolated intensities at the calculated 2theta values
|
||||
int.interp <- matrix(NA, dim(thth.a)[1], dim(thth.a)[2])
|
||||
for (j in 2:dim(int.interp)[2]) {
|
||||
int.interp[,1] <- intensity
|
||||
int.interp[,j] <- approx(thth, intensity, xout = thth.a[,j])$y
|
||||
}
|
||||
# So far, we have just replaced the old thth-scale with a new one,
|
||||
# and calculated the intensitites by linearly interpolating from the old intensities.
|
||||
|
||||
# Intensities times lever weights, P(j) (this is what you will substract from int.interp)
|
||||
int.p <- matrix(NA, length(theta), length(weights[[lever]]))
|
||||
for (j in 1:length(lengths[[lever]])) {
|
||||
int.p[,j] <- weights[[lever]][j] * int.interp[,j + 1]
|
||||
}
|
||||
|
||||
# Calculate intensities with Ka2 contribution stripped
|
||||
#int.stripped <- c(int.interp[,-1]) - rowSums(int.p)
|
||||
int.stripped <- thth - rowSums(int.p)
|
||||
#corr.tmp.df <- data.frame(thth = c(thth.a[,-1]), int.corr = int.stripped)
|
||||
corr.tmp.df <- data.frame(thth = thth, int.corr = int.stripped)
|
||||
corr.df <- corr.tmp.df[order(corr.tmp.df$thth), ]
|
||||
#row.names(corr.df) <- seq(1, length(corr.df$thth))
|
||||
|
||||
# Make a dataframe of thth.a and int.a and order it by thth
|
||||
#df.a <- data.frame(thth = c(thth.a), intensity = c(int.a))
|
||||
#df.as <- df.a[order(df.a$thth), ]
|
||||
#row.names(df.as) <- seq(1,length(df.as$thth)) # fixes row names order
|
||||
# df.as is exactly as the original data, just with the correct number of
|
||||
# interpolated thth and intensity values included.
|
||||
|
||||
|
||||
return(corr.df)
|
||||
|
||||
|
||||
|
||||
|
||||
# Perhaps thth.a and int.p.terms are the new x and y
|
||||
|
||||
|
||||
# Intensities in summation term (this is the Ka2 correction)
|
||||
#int.p <- rowSums(int.p.terms)
|
||||
|
||||
# Collapse the matrix thth.a into a single column
|
||||
#thth.ai <- matrix(NA, dim(thth.a)[1] * dim(thth.a)[2], 2)
|
||||
#thth.ai[,1] <- sort(c(thth.a))
|
||||
#thth.ai[,2] <- approx(thth, intensity, xout = thth.ai[,1])$y
|
||||
|
||||
|
||||
# Built-in functions in R to interpolate or extrapolate
|
||||
# stats::approx Linear interpolation
|
||||
# Hmisc::approxExtrap Linear extrapolation
|
||||
# go with linear functions for now, see how that works out
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# This is NOT necessary
|
||||
# This vector helps convert from lever to actual number of parameter pairs
|
||||
#parpairs <- numeric()
|
||||
#for (j in 1:length(weights)) {
|
||||
# parpairs[j] <- length(weights[[j]])
|
||||
#}
|
||||
|
||||
|
||||
##### STILL UNDER CONSTRUCTION ####
|
||||
}
|
||||
@ -1,15 +0,0 @@
|
||||
##################################################
|
||||
################ pearson.beta ####################
|
||||
##################################################
|
||||
pearson.beta <- function(m, w) {
|
||||
# Function for calculating integral breadth, \beta, of
|
||||
# an X-ray reflection using the half-width of FWHM (w)
|
||||
# and m parameter from a Pearson fit.
|
||||
## Args: m - shape parameter from Pearson fit (numeric)
|
||||
## w - half the FWHM for the reflection (numeric)
|
||||
## **** NOTE: gamma(x) out of range for x > 171
|
||||
## Values: The integral breadth of the reflection (numeric)
|
||||
# Ref: Birkholz2006, page 96 (Eq. 3.5)
|
||||
beta <- ((pi * 2^(2 * (1 - m)) * gamma(2 * m - 1)) /
|
||||
((2^(1 / m) - 1) * (gamma(m))^2)) * w
|
||||
}
|
||||
@ -1,8 +0,0 @@
|
||||
print.xtable.booktabs <- function(x){
|
||||
# Make xtable print booktabs tables
|
||||
# ARGS: x (matrix)
|
||||
require(xtable)
|
||||
print(xtable(x), floating=F, hline.after=NULL,
|
||||
add.to.row=list(pos=list(-1,0, nrow(x)),
|
||||
command=c("\\toprule ", "\\midrule ", "\\bottomrule ")))
|
||||
}
|
||||
@ -1,50 +0,0 @@
|
||||
split.muxd <- function(uxdfile, range.descriptor) {
|
||||
# Fix this some other day!!
|
||||
# Function that reads an UXD file which contains several ranges
|
||||
# (created in a programmed run, for example) and splits it into
|
||||
# its ranges and writes to that number of files
|
||||
# Arguments
|
||||
# :: uxdfile (filename with extension)
|
||||
# :: range.descriptor (an array with as many elements as
|
||||
# there are ranges in the uxdfile)
|
||||
# Returns: nothing. Writes files to HDD.
|
||||
|
||||
cchar <- "[;_]" #regexpr matching the comment characters used in Bruker's UXD
|
||||
cdata <- "[0-9]" #regexpr matching one character of any digit
|
||||
# Create filenames for the output # no longer used, return dataframe instead
|
||||
datafile <- paste(uxdfile,"-",range.descriptor,".data",sep="")
|
||||
|
||||
# Read the input multirange file
|
||||
f <- readLines(uxdfile, n=-1)
|
||||
|
||||
|
||||
# This way we identify data rows by looking for numeric characters.
|
||||
#wh <- regexpr("[0-9]", f)
|
||||
# This way we identify header rows
|
||||
# Later we will assume that all other rows are data
|
||||
wh <- regexpr(cchar, f)
|
||||
|
||||
mh <- wh[1:length(wh)] # this gives you the corresponding index vector
|
||||
# the value of each element corresponds to the position of the regexp match.
|
||||
# value = 1 means the first character of the row is cchar (row is header)
|
||||
# value =-1 means no cchar occur on the row (row is data)
|
||||
|
||||
#length(mh[mh == -1]) #total number of datarows in uxdfile
|
||||
#mh[mh > 1 | mh < 0] <- 0 #set all header-rows to zero (just to make things easier)
|
||||
|
||||
i <- seq(1, length(mh) - 1, 1)
|
||||
j <- seq(2, length(mh), 1)
|
||||
starts <- which(mh[i] == 1 & mh[j] != 1) + 1 #start indices
|
||||
ends <- which(mh[i] != 1 & mh[j] == 1) #end indices, except the last
|
||||
ends <- c(ends, length(mh)) #fixed the last index of ends
|
||||
|
||||
#ff <- data.frame(NULL)
|
||||
for (s in 1:length(range.descriptor)) {
|
||||
matrix(scan(file=textConnection(f[starts[s]:ends[s]]),
|
||||
what = numeric()), ncol=2, byrow=T)
|
||||
}
|
||||
names(ff) <- c("sampleid", "angle", "intensity")
|
||||
|
||||
# Return dataframe
|
||||
ff
|
||||
}
|
||||
@ -1,17 +0,0 @@
|
||||
##################################################
|
||||
################## strip.ka2 #####################
|
||||
##################################################
|
||||
strip.ka <- function(angle, intensity) {
|
||||
# Function that strips (removes) the Cu Ka2 contribution
|
||||
# from an experimental diffractogram (2Th, I(2Th))
|
||||
# using the Rachinger algorithm.
|
||||
## Args: 2theta vector
|
||||
## intensity vector
|
||||
## Values: intensity vector
|
||||
# Ref: Birkholz2006, page 99 (Eq. 3.12), and page 31
|
||||
Cu.Kai <- 0.154059 # nanometre
|
||||
Cu.Kaii <- 0.154441 # nanometre
|
||||
alpha.int.ratio <- 0.5
|
||||
int.stripped <- intensity - alpha.int.ratio*intensity
|
||||
return(int.stripped)
|
||||
}
|
||||
@ -1,49 +0,0 @@
|
||||
##################################################
|
||||
################### uxd2mtx ######################
|
||||
##################################################
|
||||
uxd2mtx <- function(uxdfile) {
|
||||
# Function for reading UXD files
|
||||
# Assumptions: data in two columns
|
||||
# Args: uxdfile (filename with extension)
|
||||
# Return value: matrix with two columns
|
||||
|
||||
cchar <- "[;_]" #regexpr matching the comment characters used in Bruker's UXD
|
||||
cdata <- "[0-9]" #regexpr matching one character of any digit
|
||||
|
||||
# A new file (datafile) containing only data will be created,
|
||||
# extension ".data" appended to uxdfile
|
||||
#datafile <- paste(uxdfile,".data",sep="")
|
||||
|
||||
ufile <- file(uxdfile, "r")
|
||||
f <- readLines(ufile, n=-1) #read _all_ lines from UXD file
|
||||
close(ufile)
|
||||
|
||||
# This way we identify data rows by looking for numeric characters.
|
||||
#wh <- regexpr("[0-9]", f)
|
||||
# This way we identify header rows
|
||||
# We assume that all other rows are data
|
||||
wh <- regexpr(cchar, f)
|
||||
|
||||
mh <- wh[1:length(wh)] # this gives you the corresponding index vector
|
||||
# the value of each element corresponds to the position of the regexp match.
|
||||
# value = 1 means the first character of the row is cchar (row is header)
|
||||
# value =-1 means no cchar occur on the row (row is data)
|
||||
|
||||
i <- seq(1, length(mh) - 1, 1)
|
||||
j <- seq(2, length(mh), 1)
|
||||
|
||||
starts <- which(mh[i] == 1 & mh[j] != 1) + 1
|
||||
ends <- length(mh)
|
||||
f <- f[starts:ends]
|
||||
|
||||
zz <- textConnection(f, "r")
|
||||
ff <- matrix(scan(zz, what = numeric()), ncol=2, byrow=T)
|
||||
close(zz)
|
||||
|
||||
#zz <- file(datafile, "w") #open connection to datafile
|
||||
#write.table(ff, file=datafile, row.names=F, sep=",")
|
||||
#close(zz)
|
||||
|
||||
# Return matrix
|
||||
ff
|
||||
}
|
||||
@ -1,93 +0,0 @@
|
||||
xrdpk <-
|
||||
function(xrd.exp, kerpk = 1, fitmaxiter = 50, gam = 1.0, scl.factor = 1.2, maxwdth=5.0) {
|
||||
|
||||
xrd.base <- baselinefit(xrd.exp, tau=2.5, gam=gam, scl.factor=scl.factor, maxwdth=maxwdth)
|
||||
|
||||
# This loop deals with the output from baselinefit()
|
||||
# It makes a "melted" dataframe in long form for each
|
||||
# separated peak for some baseline parameters
|
||||
xrd.pks <- data.frame()
|
||||
xrd.pks.basl <- data.frame()
|
||||
xrd.pks.pmg <- data.frame()
|
||||
xrd.pks.spl <- data.frame()
|
||||
peaks <- 1:length(xrd.base$npks)
|
||||
for (s in peaks) {
|
||||
# recorded data in long form by separated peak
|
||||
xrd.pks <- rbind(xrd.pks, # column names assigned after loop
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
xrd.exp[xrd.base$indlsep[s]:xrd.base$indrsep[s], ]))
|
||||
# the calculated baseline in long form by separated peak
|
||||
xrd.pks.basl <- rbind(xrd.pks.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = xrd.exp[xrd.base$indlsep[s]:xrd.base$indrsep[s], ],
|
||||
y = xrd.base$baseline$basisl[xrd.base$indlsep[s]:xrd.base$indrsep[s]]))
|
||||
# the taut string estimation in long form by separated peak
|
||||
xrd.pks.pmg <- rbind(xrd.pks.pmg,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = xrd.exp[xrd.base$indlsep[s]:xrd.base$indrsep[s], ],
|
||||
y = xrd.base$pmg$fn[xrd.base$indlsep[s]:xrd.base$indrsep[s]]))
|
||||
# the weighted smoothed spline in long form by separated peak
|
||||
xrd.pks.spl <- rbind(xrd.pks.spl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = xrd.exp[xrd.base$indlsep[s]:xrd.base$indrsep[s], ],
|
||||
y = xrd.base$spl$reg[xrd.base$indlsep[s]:xrd.base$indrsep[s]]))
|
||||
}
|
||||
# Column names assigned to d.pks
|
||||
names(xrd.pks) <- c("peak", "kernel", "x", "y")
|
||||
|
||||
|
||||
# This loop calls pkdecompint() on each peak separately
|
||||
# It makes a "melted" dataframe in long form for:
|
||||
xrd.fit <- list() # holds pkdecompint output
|
||||
xrd.fit.fitpk <- data.frame() # contains fitting curves
|
||||
xrd.fit.parpk <- data.frame() # physical parameters by peak and kernel
|
||||
xrd.nobasl <- data.frame() # data with baseline removed
|
||||
peaks <- 1:length(xrd.base$npks)
|
||||
for (s in peaks) {
|
||||
######## PKDECOMPINT ########
|
||||
if (length(kerpk) > 1) {
|
||||
# set number of kernels per peak manually
|
||||
xrd.fit[[s]] <- pkdecompint(xrd.base, intnum = s,
|
||||
k = kerpk[s], maxiter = fitmaxiter)
|
||||
} else {
|
||||
# use number of kernels determined by baselinefit()
|
||||
xrd.fit[[s]] <- pkdecompint(xrd.base, intnum = s,
|
||||
k = xrd.base$npks[s], maxiter = fitmaxiter)
|
||||
}
|
||||
# Setup the dataframe that makes up the peak table
|
||||
for (kernel in 1:xrd.fit[[s]]$num.ker) {
|
||||
xrd.fit.parpk <- rbind(xrd.fit.parpk,
|
||||
data.frame(peak = factor(xrd.fit[[s]]$intnr),
|
||||
kernel = factor(kernel),
|
||||
x = xrd.fit[[s]]$parpks[kernel, "loc"],
|
||||
height = xrd.fit[[s]]$parpks[kernel, "height"],
|
||||
area = xrd.fit[[s]]$parpks[kernel, "intens"],
|
||||
fwhm = xrd.fit[[s]]$parpks[kernel, "FWHM"],
|
||||
beta = xrd.fit[[s]]$parpks[kernel, "intens"] / xrd.fit[[s]]$parpks[kernel, "height"],
|
||||
m = xrd.fit[[s]]$parpks[kernel, "m"],
|
||||
accept = xrd.fit[[s]]$accept))
|
||||
xrd.fit.fitpk <- rbind(xrd.fit.fitpk,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = factor(kernel),
|
||||
x = xrd.fit[[s]]$x,
|
||||
y = xrd.fit[[s]]$fitpk[kernel, ]))
|
||||
}
|
||||
xrd.nobasl <- rbind(xrd.nobasl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
x = xrd.fit[[s]]$x,
|
||||
y = xrd.fit[[s]]$y))
|
||||
}
|
||||
|
||||
|
||||
|
||||
return(list(xrd.base = xrd.base,
|
||||
xrd.peaks = xrd.pks,
|
||||
xrd.fit.parpk = xrd.fit.parpk,
|
||||
xrd.fit.fitpk = xrd.fit.fitpk,
|
||||
xrd.nobasl = xrd.nobasl))
|
||||
|
||||
}
|
||||
@ -1,80 +0,0 @@
|
||||
xrdpkWrapper <-
|
||||
function(data.exp, run, override = FALSE,
|
||||
kerpk = 1, fitmaxiter = 50, gam = 1.0, scl.factor = 1.2, maxwdth=5.0) {
|
||||
|
||||
print("... Started xrdpkWrapper")
|
||||
|
||||
# check if xrdpk has already completed successfully for the current job
|
||||
current.dirname <- getwd()
|
||||
print(current.dirname)
|
||||
current.filename <- "xrd-peak-data.rda"
|
||||
xrddatafile <- paste(current.dirname, current.filename, sep = "/")
|
||||
|
||||
|
||||
if (file.exists(xrddatafile) && !override) {
|
||||
print("... Started if-clause 1")
|
||||
|
||||
# File already exists
|
||||
# return the data using load() or data()
|
||||
|
||||
load(file = xrddatafile)
|
||||
|
||||
if (run > length(xrdres)) {
|
||||
|
||||
print("... Started if-clause 1:1")
|
||||
|
||||
# then it does not really exist
|
||||
xrdres[[run]] <- xrdpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
save(xrdres, file = xrddatafile)
|
||||
|
||||
print("... Ended if-clause 1:1")
|
||||
}
|
||||
|
||||
print("... Ended if-clause 1")
|
||||
|
||||
return(xrdres)
|
||||
} else {
|
||||
# File does not exist
|
||||
# OR override is TRUE
|
||||
|
||||
print("... Started else-clause 1")
|
||||
|
||||
# If file does not exist at all, run all necessary code to re-create it
|
||||
if (!file.exists(xrddatafile)) {
|
||||
xrdres <- list()
|
||||
print("... xrdres list created")
|
||||
|
||||
xrdres[[run]] <-
|
||||
xrdpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
|
||||
save(xrdres, file = xrddatafile)
|
||||
} else {
|
||||
# File already exists, but override is TRUE
|
||||
load(file = xrddatafile)
|
||||
|
||||
xrdres[[run]] <-
|
||||
xrdpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
|
||||
save(xrdres, file = xrddatafile)
|
||||
}
|
||||
|
||||
print("... Ended else-clause 1")
|
||||
|
||||
return(xrdres)
|
||||
}
|
||||
}
|
||||
@ -1,104 +0,0 @@
|
||||
xrfpk <-
|
||||
function(data.exp, kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1, maxwdth = 10) {
|
||||
|
||||
print("... Starting baseline fitting")
|
||||
|
||||
data.basl <- baselinefit(data.exp,
|
||||
tau = 2.0,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
|
||||
print("... Ended baseline fitting")
|
||||
|
||||
# This loop deals with the output from baselinefit()
|
||||
# It makes a "melted" dataframe in long form for each
|
||||
# separated peak for some baseline parameters
|
||||
data.pks <- data.frame()
|
||||
data.pks.basl <- data.frame()
|
||||
data.pks.pmg <- data.frame()
|
||||
data.pks.spl <- data.frame()
|
||||
peaks <- 1:length(data.basl$npks)
|
||||
|
||||
|
||||
for (s in peaks) {
|
||||
# recorded data in long form by separated peak
|
||||
data.pks <- rbind(data.pks, # column names assigned after loop
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], ]))
|
||||
# the calculated baseline in long form by separated peak
|
||||
data.pks.basl <- rbind(data.pks.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$baseline$basisl[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
# the taut string estimation in long form by separated peak
|
||||
data.pks.pmg <- rbind(data.pks.pmg,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$pmg$fn[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
# the weighted smoothed spline in long form by separated peak
|
||||
data.pks.spl <- rbind(data.pks.spl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = NA,
|
||||
x = data.exp[data.basl$indlsep[s]:data.basl$indrsep[s], 1],
|
||||
y = data.basl$spl$reg[data.basl$indlsep[s]:data.basl$indrsep[s]]))
|
||||
}
|
||||
|
||||
|
||||
|
||||
# Column names assigned to d.pks
|
||||
names(data.pks) <- c("peak", "kernel", "x", "y")
|
||||
|
||||
|
||||
# This loop calls pkdecompint() on each peak separately
|
||||
# It makes a "melted" dataframe in long form for:
|
||||
data.fit <- list() # holds pkdecompint output
|
||||
data.fit.fitpk <- data.frame() # contains fitting curves
|
||||
data.fit.parpk <- data.frame() # physical parameters by peak and kernel
|
||||
data.fit.basl <- data.frame() # data with baseline removed
|
||||
peaks <- 1:length(data.basl$npks)
|
||||
for (s in peaks) {
|
||||
######## PKDECOMPINT ########
|
||||
if (length(kerpk) > 1) {
|
||||
# set number of kernels per peak manually
|
||||
data.fit[[s]] <- pkdecompint(data.basl, intnum = s,
|
||||
k = kerpk[s], maxiter = fitmaxiter)
|
||||
} else {
|
||||
# use number of kernels determined by baselinefit()
|
||||
data.fit[[s]] <- pkdecompint(data.basl, intnum = s,
|
||||
k = data.basl$npks[s], maxiter = fitmaxiter)
|
||||
}
|
||||
# Setup the dataframe that makes up the peak table
|
||||
for (kernel in 1:data.fit[[s]]$num.ker) {
|
||||
data.fit.parpk <- rbind(data.fit.parpk,
|
||||
data.frame(peak = factor(data.fit[[s]]$intnr),
|
||||
kernel = factor(kernel),
|
||||
x = data.fit[[s]]$parpks[kernel, "loc"],
|
||||
height = data.fit[[s]]$parpks[kernel, "height"],
|
||||
area = data.fit[[s]]$parpks[kernel, "intens"],
|
||||
fwhm = data.fit[[s]]$parpks[kernel, "FWHM"],
|
||||
m = data.fit[[s]]$parpks[kernel, "m"],
|
||||
accept = data.fit[[s]]$accept))
|
||||
data.fit.fitpk <- rbind(data.fit.fitpk,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
kernel = factor(kernel),
|
||||
x = data.fit[[s]]$x,
|
||||
y = data.fit[[s]]$fitpk[kernel, ]))
|
||||
}
|
||||
data.fit.basl <- rbind(data.fit.basl,
|
||||
data.frame(peak = factor(peaks[s]),
|
||||
x = data.fit[[s]]$x,
|
||||
y = data.fit[[s]]$y))
|
||||
}
|
||||
|
||||
|
||||
|
||||
return(list(data.basl = data.basl,
|
||||
data.peaks = data.pks,
|
||||
data.fit.parpk = data.fit.parpk,
|
||||
data.fit.fitpk = data.fit.fitpk,
|
||||
data.fit.basl = data.fit.basl))
|
||||
}
|
||||
@ -1,67 +0,0 @@
|
||||
xrfpkWrapper <-
|
||||
function(data.exp, run, override = FALSE,
|
||||
kerpk = 1, fitmaxiter = 100, gam = 0.64, scl.factor = 0.06, maxwdth = 10) {
|
||||
# kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor = 0.1) {
|
||||
# the override flag is IN USE
|
||||
|
||||
print("... Started xrfpkWrapper")
|
||||
|
||||
# check if xrfpk has already completed successfully for the current job
|
||||
current.dirname <- getwd()
|
||||
print(current.dirname)
|
||||
current.filename <- "xrf-peak-data.rda"
|
||||
xrfdatafile <- paste(current.dirname, current.filename, sep = "/")
|
||||
|
||||
|
||||
|
||||
if (file.exists(xrfdatafile) && !override) {
|
||||
print("... Started if-clause 1")
|
||||
|
||||
# File already exists
|
||||
# return the data using load() or data()
|
||||
|
||||
load(file = xrfdatafile)
|
||||
|
||||
if (run > length(xrfres)) {
|
||||
|
||||
print("... Started if-clause 1:1")
|
||||
|
||||
# the it does not really exist
|
||||
xrfres[[run]] <- xrfpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
save(xrfres, file = xrfdatafile)
|
||||
|
||||
print("... Ended if-clause 1:1")
|
||||
}
|
||||
|
||||
print("... Ended if-clause 1")
|
||||
|
||||
return(xrfres)
|
||||
} else {
|
||||
|
||||
print("... Started else-clause 1")
|
||||
|
||||
if (!exists("xrfres")) {
|
||||
xrfres <- list()
|
||||
print("... xrfres list created")
|
||||
}
|
||||
|
||||
# Need to call xrfpk() and save its results to file as above
|
||||
xrfres[[run]] <- xrfpk(data.exp,
|
||||
kerpk = kerpk,
|
||||
fitmaxiter = fitmaxiter,
|
||||
gam = gam,
|
||||
scl.factor = scl.factor,
|
||||
maxwdth = maxwdth)
|
||||
save(xrfres, file = xrfdatafile)
|
||||
|
||||
print("... Ended else-clause 1")
|
||||
|
||||
return(xrfres)
|
||||
}
|
||||
|
||||
}
|
||||
@ -1,3 +0,0 @@
|
||||
I have been attempting to fine-tune the parameters of the \texttt{baselinefit()} function.
|
||||
|
||||
Increasing the number of iterations above 100 have been useless. Small changes in \texttt{gam} leads to noticeable changes in the peak and kernel distribution.
|
||||
@ -1,205 +0,0 @@
|
||||
source("/home/taha/chepec/chetex/common/R/common/ProvideSampleId.R")
|
||||
|
||||
xrfspectro2df <- function(smpfile) {
|
||||
## Description:
|
||||
## Total remake of xrfspectro2df(). Idea is to accomodate all 6 possible
|
||||
## datasets of each SMP file, plus the attributes.
|
||||
## Reads XRF textfile from XLAB SPECTRO XRF.
|
||||
## Stores data in data frame and parameters in an attributed dataframe.
|
||||
## Usage:
|
||||
## xrfspectro2df(smpfile)
|
||||
## Arguments:
|
||||
## smpfile: character string, the full filename
|
||||
## (with path) to one SMP file (ASCII).
|
||||
## Value:
|
||||
## A dataframe with attributed dataframe
|
||||
|
||||
#### ONLY BOTHER WITH THE FIRST MEASUREMENT IN THE SMP-FILE.
|
||||
|
||||
filecon <- file(smpfile, "r")
|
||||
smpcontents <- readLines(filecon, n = -1) #read all lines of input file
|
||||
close(filecon)
|
||||
|
||||
# Parameter table
|
||||
# Those are the parameter we may access later in this function
|
||||
xrf.param <- data.frame(stringsAsFactors = FALSE,
|
||||
matrix(c("Method", "^Method:",
|
||||
"Job", "^Job:",
|
||||
"Status", "^Status:",
|
||||
"Description", "^Description:",
|
||||
"Date", "^Date\\sof\\sMeasurement:",
|
||||
"Measurements", "^Measurements:",
|
||||
"Voltage", "^Voltage:",
|
||||
"Current", "^Current:",
|
||||
"Target", "^Target:",
|
||||
"Duration", "^Meas\\.\\sDuration:",
|
||||
"Impulse", "^Imp\\.\\sRate:",
|
||||
"DeadTime", "^Rel\\.\\sDead\\sTime:",
|
||||
"FirstChannel", "^First\\sChannel:",
|
||||
"LastChannel", "^Last\\sChannel:",
|
||||
"PeakTime", "^Peak\\sTime:",
|
||||
"Gain", "^Gain:",
|
||||
"ZeroPeak", "^Zero\\sPeak:",
|
||||
"Data", "^Kanal\\s[\\d]+:"),
|
||||
ncol = 2, byrow = T))
|
||||
names(xrf.param) <- c("parameter", "regexp")
|
||||
|
||||
|
||||
# Data table
|
||||
# Contains the regexp used for identifiying rows containing data
|
||||
xrf.data <- data.frame(stringsAsFactors = FALSE,
|
||||
matrix(c("Data", "^Kanal\\s[\\d]+:"), ncol = 2, byrow = T))
|
||||
names(xrf.data) <- c("parameter", "regexp")
|
||||
|
||||
|
||||
|
||||
# Find out how many measurements we have in this
|
||||
# file by accessing the Measurements field
|
||||
n_measurements <- as.numeric(strsplit(gsub("^\\t", "",
|
||||
strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Measurements", select = "regexp")$regexp,
|
||||
smpcontents) == 1)], ":")[[1]][2]), "\\t")[[1]])
|
||||
# If more than one measurement, issue warning
|
||||
if (n_measurements > 1) {
|
||||
warning(paste(paste(n_measurements, " measurements detected in ",
|
||||
basename(smpfile), sep = ""),
|
||||
"Only the first measurement will be recorded",
|
||||
sep = "\n", collapse = ""))
|
||||
}
|
||||
|
||||
|
||||
# How many rows of data?
|
||||
n_rowsdata <-
|
||||
length(which(regexpr(subset(xrf.data, parameter == "Data",
|
||||
select = "regexp")$regexp, smpcontents, perl = TRUE) == 1))
|
||||
|
||||
# Build an empty matrix big enough to hold all data
|
||||
# (i.e., ncol = 3, and nrow = n_rowsdata * n_measurements)
|
||||
data.long <- data.frame(matrix(NA, ncol = 5, nrow = 6 * n_rowsdata))
|
||||
names(data.long) <- c("sampleid", "measurement", "channel", "energy", "counts")
|
||||
|
||||
|
||||
data.mtx <- matrix(NA, ncol = 6, nrow = n_rowsdata)
|
||||
for (j in 1:n_rowsdata) {
|
||||
data.mtx[j, ] <- as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.data, parameter == "Data",
|
||||
select = "regexp")$regexp, smpcontents, perl = TRUE) == 1)], ":\\t")[[j]][2], "\\t")[[1]])
|
||||
}
|
||||
|
||||
|
||||
# Sampleid to column 1
|
||||
data.long[, 1] <- rep(ProvideSampleId(smpfile), n_rowsdata)
|
||||
# Channel to column 3
|
||||
data.long[, 3] <- rep(seq(1, n_rowsdata), dim(data.mtx)[2])
|
||||
for (c in 1:6) {
|
||||
# Measurement no. in column 2
|
||||
data.long[((c * n_rowsdata) - n_rowsdata + 1):(((c + 1) * n_rowsdata) - n_rowsdata), 2] <- rep(c, n_rowsdata)
|
||||
# Counts in column 5
|
||||
data.long[((c * n_rowsdata) - n_rowsdata + 1):(((c + 1) * n_rowsdata) - n_rowsdata), 5] <- data.mtx[, c]
|
||||
}
|
||||
|
||||
# Drop all rows with measurement-number not equal to 1
|
||||
data.long <- subset(data.long, measurement == 1)
|
||||
|
||||
|
||||
# Fetch the measurement parameters
|
||||
data.long[, subset(xrf.param, parameter == "Date")$parameter] <-
|
||||
rep(substr(strsplit(gsub("^\\t", "", strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Date")$regexp, smpcontents) == 1)], ":")[[1]][2]),
|
||||
"\\t")[[1]][1], 1, 8), n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Method")$parameter] <-
|
||||
rep(strsplit(gsub("^\\t", "", strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Method")$regexp, smpcontents) == 1)], ":")[[1]][2]),
|
||||
"\\t")[[1]][1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Job")$parameter] <-
|
||||
rep(strsplit(gsub("^\\t", "", strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Job")$regexp, smpcontents) == 1)], ":")[[1]][2]),
|
||||
"\\t")[[1]][1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Status")$parameter] <-
|
||||
rep(strsplit(gsub("^\\t", "", strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Status")$regexp, smpcontents) == 1)], ":")[[1]][2]),
|
||||
"\\t")[[1]][1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Description")$parameter] <-
|
||||
rep(gsub("^\\t", "", strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Description")$regexp, smpcontents) == 1)], ":")[[1]][2]),
|
||||
n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Voltage")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Voltage")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Current")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Current")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Target")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Target")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Duration")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Duration")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Impulse")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Impulse")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "DeadTime")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "DeadTime")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "FirstChannel")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "FirstChannel")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "LastChannel")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "LastChannel")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "PeakTime")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "PeakTime")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "Gain")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "Gain")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
data.long[, subset(xrf.param, parameter == "ZeroPeak")$parameter] <-
|
||||
rep(as.numeric(strsplit(strsplit(smpcontents[which(regexpr(subset(xrf.param,
|
||||
parameter == "ZeroPeak")$regexp, smpcontents, perl = TRUE) == 1)],
|
||||
":\\t")[[1]][2], "\\t")[[1]])[1], n_rowsdata)
|
||||
|
||||
|
||||
# Convert channel into energy scale
|
||||
# Using the following assumptions:
|
||||
# 1. Zero peak is always the strongest (highest) peak in the spectrum
|
||||
# The channel with maximum counts should correspond to 0 keV
|
||||
# This gives a one-channel deviation from what the instrument shows
|
||||
# for a 12.5 keV range measurement using 1024 channels (so far)
|
||||
# This is good enough for our purposes, since the peak energies for most
|
||||
# ions do not match with reference values without a correction term anyway.
|
||||
max.channel <- which(data.long$counts == max(data.long$counts))
|
||||
data.long$energy <- (data.long$channel * (data.long$Gain / data.long$LastChannel)) -
|
||||
((max.channel / data.long$LastChannel) * data.long$Gain)
|
||||
# Save the maxchannel to the returned dataframe
|
||||
data.long$ZeroChannel <- rep(max.channel, n_rowsdata)
|
||||
|
||||
# Calculate energy from channel # this is no longer viable
|
||||
#data.long$energy <- (data.long$channel * (data.long$Gain / data.long$LastChannel)) -
|
||||
# ((24 / data.long$LastChannel) * data.long$Gain)
|
||||
|
||||
return(data.long)
|
||||
}
|
||||
@ -1,19 +0,0 @@
|
||||
# common.R
|
||||
# General-purpose functions
|
||||
# Taha Ahmed, Jan 2011
|
||||
|
||||
# CONTENTS Status Depends on
|
||||
# -------- ------ ----------
|
||||
# >>>> LinearBaseline deprecated ?
|
||||
# >>>> int2padstr ?
|
||||
# >>>> It2charge deprecated ?
|
||||
# >>>> ProvideSampleId ?
|
||||
# >>>> Celsius2Kelvin ?
|
||||
# >>>> Kelvin2Celsius ?
|
||||
# >>>> as.radians ?
|
||||
# >>>> as.degrees ?
|
||||
# >>>> molarity2mass ?
|
||||
# >>>> AVS2SHE -
|
||||
# >>>> SHE2AVS -
|
||||
# >>>> ConvertRefPotEC -
|
||||
# >>>> ConvertRefPot ConvertRefPotEC, SHE2AVS, AVS2SHE
|
||||
@ -1,35 +0,0 @@
|
||||
##################################################
|
||||
################### It2charge ####################
|
||||
##################################################
|
||||
It2charge <- function (time, current) {
|
||||
## **** STOP USING THIS FUNCTION *** CAUSED WEIRD, UNREPRODUCIBLE ERRORS /110304
|
||||
## Description:
|
||||
## Calculates cumulative charge, differentials, etc. from
|
||||
## amperometric data (current and time).
|
||||
## __Intended to be used from within other functions (CHI.R)__
|
||||
## Usage:
|
||||
## It2charge(time, current)
|
||||
## Arguments:
|
||||
## time: a vector with time data.
|
||||
## current: a vector of the same length as time, with current data.
|
||||
## May be either currents or current densities, no matter.
|
||||
## Value:
|
||||
## Returns a dataframe with columns:
|
||||
## timediff, dIdt, charge, sumcharge
|
||||
#
|
||||
# Calculate the time vector difference
|
||||
timediff <- c(time[1], diff(time))
|
||||
# timediff times the current gives the charge,
|
||||
# since the time vector can be considered as
|
||||
# the cumulative time, while we need to multiply
|
||||
# the current with the time elapsed since the last
|
||||
# current measurement (the timediff).
|
||||
charge <- current * timediff
|
||||
dIdt <- current / time
|
||||
# Return value
|
||||
ff <- data.frame(timediff = timediff,
|
||||
dIdt = dIdt, charge = charge,
|
||||
# perhaps it is more correct to calculate cumsum of the absolute charge?
|
||||
sumcharge = cumsum(charge))
|
||||
return(ff)
|
||||
}
|
||||
@ -1,23 +0,0 @@
|
||||
##################################################
|
||||
################ LinearBaseline ##################
|
||||
##################################################
|
||||
LinearBaseline <- function (potential, current, iplim) {
|
||||
## Arguments:
|
||||
## potential: full half-cycle, potentials
|
||||
## current: full half-cycle, currents
|
||||
## iplim: interpolation limits along x (potential)
|
||||
## Value:
|
||||
## A dataframe with two columns, potential and current
|
||||
## (of the calculated baseline)
|
||||
# Construct potential-current dataframe
|
||||
sweep <- data.frame(potential = potential, current = current)
|
||||
#
|
||||
sweep.iplim <- subset(subset(sweep, potential > iplim[1]), potential < iplim[2])
|
||||
sweep.baseline <- data.frame(potential = approxExtrap(sweep.iplim$potential,
|
||||
sweep.iplim$current, xout = sweep$potential, method = "linear")$x,
|
||||
current = approxExtrap(sweep.iplim$potential, sweep.iplim$current,
|
||||
xout = sweep$potential, method = "linear")$y)
|
||||
sweep.data <- data.frame(potential = sweep.baseline$potential,
|
||||
current = sweep.baseline$current)
|
||||
return(sweep.data)
|
||||
}
|
||||
@ -1,10 +0,0 @@
|
||||
# To source a bunch of files in the same directory
|
||||
|
||||
sourceDir <- function(path, trace = TRUE) {
|
||||
lsDir <- list.files(path, pattern = "\\.[Rr]$")
|
||||
for (i in lsDir) {
|
||||
if(trace) {cat(i, ":")}
|
||||
source(file.path(path, i))
|
||||
if(trace) {cat("\n")}
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue