diff --git a/AUTOLAB/AUTOLAB.info b/AUTOLAB/AUTOLAB.info new file mode 100644 index 0000000..5332a72 --- /dev/null +++ b/AUTOLAB/AUTOLAB.info @@ -0,0 +1,7 @@ +# AUTOLAB.R +# Functions to read and manipulate data from AUTOLAB potentiostats/galvanostats +# Taha Ahmed, June 2011 + +# CONTENTS +# >>>> amp2df +# >>>> ocp2df diff --git a/AUTOLAB.R b/AUTOLAB/amp2df.R similarity index 58% rename from AUTOLAB.R rename to AUTOLAB/amp2df.R index 6da256b..0eeaac9 100644 --- a/AUTOLAB.R +++ b/AUTOLAB/amp2df.R @@ -1,14 +1,4 @@ -# AUTOLAB.R -# Functions to read and manipulate data from AUTOLAB potentiostats/galvanostats -# Taha Ahmed, June 2011 - -# CONTENTS source("/home/taha/chepec/chetex/common/R/common.R") -# >>>> amp2df -# >>>> ocp2df - - - ################################################## #################### amp2df ###################### @@ -96,67 +86,4 @@ amp2df <- function(datafilename, wearea = 1) { chargedensity = chargedensity) # return(ff) -} - - - - -################################################## -#################### 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) } \ No newline at end of file diff --git a/AUTOLAB/ocp2df.R b/AUTOLAB/ocp2df.R new file mode 100644 index 0000000..38e6dae --- /dev/null +++ b/AUTOLAB/ocp2df.R @@ -0,0 +1,61 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI.R b/CHI.R deleted file mode 100644 index 1be8745..0000000 --- a/CHI.R +++ /dev/null @@ -1,785 +0,0 @@ -# CHI.R -# Functions to read and manipulate data from the CHI760 potentiostat/galvanostat -# Taha Ahmed, Jan 2011 - Feb 2011 - -# CONTENTS -source("/home/taha/chepec/chetex/common/R/common.R") -# >>>> mps2df -# >>>> ocp2df -# >>>> chronocm2df -# >>>> chronoamp2df -# >>>> amperometry2df -# >>>> cv2df -# >>>> lsv2df - - - -################################################## -################### 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) -} - - - - -################################################## -################### 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) -} - - - -################################################## -################# 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) -} - - - -################################################## -################# 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) -} - - - -################################################## -############### 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) -} - - - -################################################## -#################### 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) -} - - - - -################################################## -################### 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) -} diff --git a/CHI/CHI.info b/CHI/CHI.info new file mode 100644 index 0000000..4d16ae0 --- /dev/null +++ b/CHI/CHI.info @@ -0,0 +1,12 @@ +# 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 diff --git a/CHI/amperometry2df.R b/CHI/amperometry2df.R new file mode 100644 index 0000000..ccee1ed --- /dev/null +++ b/CHI/amperometry2df.R @@ -0,0 +1,123 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/chronoamp2df.R b/CHI/chronoamp2df.R new file mode 100644 index 0000000..f3157b9 --- /dev/null +++ b/CHI/chronoamp2df.R @@ -0,0 +1,122 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/chronocm2df.R b/CHI/chronocm2df.R new file mode 100644 index 0000000..12550dd --- /dev/null +++ b/CHI/chronocm2df.R @@ -0,0 +1,78 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/cv2df.R b/CHI/cv2df.R new file mode 100644 index 0000000..c209996 --- /dev/null +++ b/CHI/cv2df.R @@ -0,0 +1,94 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/lsv2df.R b/CHI/lsv2df.R new file mode 100644 index 0000000..0facc0a --- /dev/null +++ b/CHI/lsv2df.R @@ -0,0 +1,114 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/mps2df.R b/CHI/mps2df.R new file mode 100644 index 0000000..ad2edd1 --- /dev/null +++ b/CHI/mps2df.R @@ -0,0 +1,147 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/CHI/ocp2df.R b/CHI/ocp2df.R new file mode 100644 index 0000000..3aef058 --- /dev/null +++ b/CHI/ocp2df.R @@ -0,0 +1,85 @@ +source("/home/taha/chepec/chetex/common/R/common.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) +} \ No newline at end of file diff --git a/INCA/INCA.info b/INCA/INCA.info new file mode 100644 index 0000000..102c59e --- /dev/null +++ b/INCA/INCA.info @@ -0,0 +1,8 @@ +# 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 diff --git a/INCA/eds2df.R b/INCA/eds2df.R new file mode 100644 index 0000000..14b1497 --- /dev/null +++ b/INCA/eds2df.R @@ -0,0 +1,88 @@ +################################################## +################### eds2df ####################### +################################################## +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) +} \ No newline at end of file diff --git a/INCA.R b/INCA/edspk.R similarity index 53% rename from INCA.R rename to INCA/edspk.R index 5eb83ea..e4b4f9c 100644 --- a/INCA.R +++ b/INCA/edspk.R @@ -1,107 +1,3 @@ -# 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 - - - - -################################################## -################### eds2df ####################### -################################################## -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) -} - - - - ################################################## #################### edspk ####################### ################################################## @@ -195,18 +91,4 @@ edspk <- function(eds.exp, kerpk = 1, fitmaxiter = 50) { eds.fit.fitpk = eds.fit.fitpk, eds.nobasl = eds.nobasl)) -} - - - - - - - - - - - - - - +} \ No newline at end of file diff --git a/LEO1550/LEO1550.info b/LEO1550/LEO1550.info new file mode 100644 index 0000000..932f5d2 --- /dev/null +++ b/LEO1550/LEO1550.info @@ -0,0 +1,6 @@ +# LEO1550.R +# Functions to read and manipulate data from the LEO1550 scanning electron microscope +# Taha Ahmed, March 2011 + +# CONTENTS +# >>>> tifftags2df diff --git a/LEO1550.R b/LEO1550/tifftags2df.R similarity index 96% rename from LEO1550.R rename to LEO1550/tifftags2df.R index d68140b..38a8131 100644 --- a/LEO1550.R +++ b/LEO1550/tifftags2df.R @@ -1,13 +1,4 @@ -# LEO1550.R -# Functions to read and manipulate data from the LEO1550 scanning electron microscope -# Taha Ahmed, March 2011 - -# CONTENTS source("/home/taha/chepec/chetex/common/R/common.R") -# >>>> tifftags2df - - - ################################################## ################ tifftags2df ##################### @@ -136,23 +127,4 @@ tifftags2df <- function(tiffimage) { } # return(tagsdf) -} - - - - - - - - - - - - - - - - - - - +} \ No newline at end of file diff --git a/Renishaw/Raman2df.R b/Renishaw/Raman2df.R new file mode 100644 index 0000000..a6678e3 --- /dev/null +++ b/Renishaw/Raman2df.R @@ -0,0 +1,34 @@ +source("/home/taha/chepec/chetex/common/R/common.R") + +################################################## +################### Raman2df ####################### +################################################## +Raman2df <- function(datafilename) { + # Function description: for reading Raman spectrum into dataframe + # + 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", "shift", "counts") + # Re-order by increasing shift + ff <- ff[order(ff$shift), ] + # And fix the row.names + row.names(ff) <- seq(1, dim(ff)[1]) + # Do not re-calculate the spectrum with evenly spaced points here! + # You must first remove cosmic peaks, and as long as that is done + # manually, re-calculation to evenly spaced shifts must also be + # done manually. + ## + return(ff) +} \ No newline at end of file diff --git a/Renishaw.R b/Renishaw/Ramanpk.R similarity index 77% rename from Renishaw.R rename to Renishaw/Ramanpk.R index 3d4deb7..5e15de6 100644 --- a/Renishaw.R +++ b/Renishaw/Ramanpk.R @@ -1,49 +1,4 @@ -# Renishaw.R -# Functions to read data from the Renishaw Raman spectrometer -# Taha Ahmed, Feb 2011 - -# CONTENTS source("/home/taha/chepec/chetex/common/R/common.R") -# >>>> Raman2df -# >>>> Ramanpk - - - - - -################################################## -################### Raman2df ####################### -################################################## -Raman2df <- function(datafilename) { - # Function description: for reading Raman spectrum into dataframe - # - 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", "shift", "counts") - # Re-order by increasing shift - ff <- ff[order(ff$shift), ] - # And fix the row.names - row.names(ff) <- seq(1, dim(ff)[1]) - # Do not re-calculate the spectrum with evenly spaced points here! - # You must first remove cosmic peaks, and as long as that is done - # manually, re-calculation to evenly spaced shifts must also be - # done manually. - ## - return(ff) -} - ################################################## ################## Ramanpk ####################### @@ -137,4 +92,4 @@ Ramanpk <- function(Raman.exp, kerpk = 1, fitmaxiter = 50, gam = 0.6, scl.factor Raman.fit.parpk = Raman.fit.parpk, Raman.fit.fitpk = Raman.fit.fitpk, Raman.nobasl = Raman.nobasl)) -} +} \ No newline at end of file diff --git a/Renishaw/Renishaw.info b/Renishaw/Renishaw.info new file mode 100644 index 0000000..38a5e32 --- /dev/null +++ b/Renishaw/Renishaw.info @@ -0,0 +1,7 @@ +# Renishaw.R +# Functions to read data from the Renishaw Raman spectrometer +# Taha Ahmed, Feb 2011 + +# CONTENTS +# >>>> Raman2df +# >>>> Ramanpk diff --git a/XRD-TF/XRD-TF.info b/XRD-TF/XRD-TF.info new file mode 100644 index 0000000..0afd3a3 --- /dev/null +++ b/XRD-TF/XRD-TF.info @@ -0,0 +1,21 @@ +# 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 diff --git a/XRD-TF/matchpdf.R b/XRD-TF/matchpdf.R new file mode 100644 index 0000000..6eb36dc --- /dev/null +++ b/XRD-TF/matchpdf.R @@ -0,0 +1,76 @@ +################################################## +################## 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) +} diff --git a/XRD-TF/muxd2df.R b/XRD-TF/muxd2df.R new file mode 100644 index 0000000..807132e --- /dev/null +++ b/XRD-TF/muxd2df.R @@ -0,0 +1,54 @@ +################################################## +################### muxd2df ###################### +################################################## +muxd2df <- function(uxdfile, range.descriptor) { + # Function that reads an UXD file which contains several ranges + # (created in a programmed run, for example) + # Arguments + # :: uxdfile (filename with extension) + # :: range.descriptor (an array with as many elements as + # there are ranges in the uxdfile) + # Returns: dataframe with 3 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) + + # 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)) { + zz <- textConnection(f[starts[s]:ends[s]], "r") + ff <- rbind(ff, data.frame(range.descriptor[s], + matrix(scan(zz, what = numeric()), ncol=2, byrow=T))) + close(zz) + } + names(ff) <- c("sampleid", "angle", "intensity") + + # Return dataframe + ff +} diff --git a/XRD-TF/muxd2ls.R b/XRD-TF/muxd2ls.R new file mode 100644 index 0000000..5419c7d --- /dev/null +++ b/XRD-TF/muxd2ls.R @@ -0,0 +1,38 @@ +################################################## +################### 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) +} diff --git a/XRD-TF/muxd2mtx.R b/XRD-TF/muxd2mtx.R new file mode 100644 index 0000000..6a8d813 --- /dev/null +++ b/XRD-TF/muxd2mtx.R @@ -0,0 +1,57 @@ +################################################## +################### 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 +} diff --git a/XRD-TF/pdf2df.R b/XRD-TF/pdf2df.R new file mode 100644 index 0000000..ab5a15c --- /dev/null +++ b/XRD-TF/pdf2df.R @@ -0,0 +1,71 @@ +################################################## +################### pdf2df ####################### +################################################## +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, this function simply strips the less-than character. + # (Perhaps not true, see the int.Tex column) + # 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) + # attr: This function sets the following attributes: + # ApplicationName, + # ApplicationVersion, + # pdfNumber, + # 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 = "") + )) + } + # + 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) +} diff --git a/XRD-TF/scherrer.R b/XRD-TF/scherrer.R new file mode 100644 index 0000000..3fb60dc --- /dev/null +++ b/XRD-TF/scherrer.R @@ -0,0 +1,15 @@ +################################################## +################## 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(), source("/home/taha/chepec/chetex/common/R/common.R") + D <- (shapeconstant * wavelength) / (as.radians(integralbreadth) * cos(as.radians(thth))) + # cos() - angles must be in radians, not degrees! + return(D) #units of angstrom +} diff --git a/XRD-TF/unborn/EliminateKa2() b/XRD-TF/unborn/EliminateKa2() new file mode 100644 index 0000000..f2c37a6 --- /dev/null +++ b/XRD-TF/unborn/EliminateKa2() @@ -0,0 +1,215 @@ +################################################## +################ 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 #### +} diff --git a/XRD-TF/unborn/pearson.beta() b/XRD-TF/unborn/pearson.beta() new file mode 100644 index 0000000..99ab4ca --- /dev/null +++ b/XRD-TF/unborn/pearson.beta() @@ -0,0 +1,15 @@ +################################################## +################ 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 +} diff --git a/XRD-TF/unborn/print.xtable.booktabs() b/XRD-TF/unborn/print.xtable.booktabs() new file mode 100644 index 0000000..02ccfe8 --- /dev/null +++ b/XRD-TF/unborn/print.xtable.booktabs() @@ -0,0 +1,8 @@ +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 "))) +} diff --git a/XRD-TF/unborn/split.muxd() b/XRD-TF/unborn/split.muxd() new file mode 100644 index 0000000..e498d2e --- /dev/null +++ b/XRD-TF/unborn/split.muxd() @@ -0,0 +1,50 @@ +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 +} diff --git a/XRD-TF/unborn/strip.ka2() b/XRD-TF/unborn/strip.ka2() new file mode 100644 index 0000000..0ebf99c --- /dev/null +++ b/XRD-TF/unborn/strip.ka2() @@ -0,0 +1,17 @@ +################################################## +################## 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) +} diff --git a/XRD-TF/uxd2df.R b/XRD-TF/uxd2df.R new file mode 100644 index 0000000..8d6fd5e --- /dev/null +++ b/XRD-TF/uxd2df.R @@ -0,0 +1,63 @@ +################################################## +#################### uxd2df ###################### +################################################## +uxd2df <- function(uxdfile) { + # Function for reading UXD files # Assumptions: data in two columns + # Args: uxdfile (filename with extension) + # Returns: dataframe with three 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] + + rgxp.sampleid <- "[^/]*(?=\\.\\w*)" ## THIS REQUIRES perl=TRUE + # Regular expression that extracts the filename out of a full path. + # Matches and extracts everything from the last forward slash (assuming Unix slashes) + # up until a dot folllowed by an arbitrary number of alphanumeric characters. + sampleidmtch <- regexpr(rgxp.sampleid, uxdfile, perl=TRUE) + # Check that there was a match + if (sampleidmtch < 0) { + # -1 means no match + sampleid <- uxdfile + # If match was unsuccessful we use the argument as passed to this function as sampleid + } + sampleid <- substr(uxdfile, sampleidmtch, (sampleidmtch + attr(sampleidmtch, "match.length") - 1)) + + zz <- textConnection(f, "r") + ff <- data.frame(sampleid, matrix(scan(zz, + what = numeric()), ncol=2, byrow=T)) + names(ff) <- c("sampleid", "angle", "intensity") + close(zz) + + #zz <- file(datafile, "w") #open connection to datafile + #write.table(ff, file=datafile, row.names=F, sep=",") + #close(zz) + + # Return dataframe + ff +} diff --git a/XRD-TF/uxd2mtx.R b/XRD-TF/uxd2mtx.R new file mode 100644 index 0000000..1256661 --- /dev/null +++ b/XRD-TF/uxd2mtx.R @@ -0,0 +1,49 @@ +################################################## +################### 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 +} diff --git a/common.R b/common.R deleted file mode 100644 index 7a94146..0000000 --- a/common.R +++ /dev/null @@ -1,353 +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 - - - -################################################## -################ 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) -} - - - - -################################################## -################## int2padstr #################### -################################################## -int2padstr <- function (ii, pchr, w) { - ## Description: - ## Converts an integer or a vector of integers to - ## a string padded with characters. - ## Usage: - ## int2padstr(ii, pchr, w) - ## Arguments: - ## ii: integer or vector of integers - ## pchr: a padding character (e.g., "0") - ## w: width of the return string (an integer) - ## Make sure to set the width longer than - ## or equal to the length of the biggest integer. - ## For example, if the integers (ii) are - ## in the range 1 - 100, set w to at least 3. - ## Value: - ## A character string or a vector of character strings - gsub(" ", pchr, formatC(ii, format="s", mode="character", width = w)) -} - - -################################################## -################### 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) -} - - - -################################################## -################ ProvideSampleId ################# -################################################## -ProvideSampleId <- function (fullpathwithfilename) { - ### OBS! Only very rudimentary error-checking. - ### If the filename is formatted as \w*-\w*-\w*, we use the middle segment, - ### otherwise we use the whole string (excluding the extension) - # Extract the name of the parent directory of the datafilename argument - substrateid <- basename(dirname(fullpathwithfilename)) - # Extract the name of the method from the filename-part - # First split the filename over all hyphens - nameparts <- strsplit(basename(fullpathwithfilename), "-")[[1]] - # If the number of nameparts exceed 3, save the whole filename as methodid, otherwise use the middle part - if (length(nameparts) > 3) { - # We need to lose the file extension from the last namepart - nameparts[length(nameparts)] <- strsplit(nameparts[length(nameparts)], "\\.")[[1]][1] - methodid <- paste(nameparts, collapse = "-") - } else { - methodid <- nameparts[2] - } - # Make an informative sampleid - sampleid <- paste(substrateid, methodid, sep = "-") - # - return(sampleid) -} - - - -################################################## -############### Celsius2Kelvin ################### -################################################## -Celsius2Kelvin <- function(Celsius) { - # Converts temperature from Celsius to Kelvin - # - # Check and correct for values below -273.15 - if (Celsius < -273.15) { - # If Celsis is less than absolute zero, set it to absolute zero - Celsius <- -273.15 - } - Kelvin <- Celsius + 273.15 - return(Kelvin) -} - - -################################################## -############### Kelvin2Celsius ################### -################################################## -Kelvin2Celsius <- function(Kelvin) { - # Converts temperature from Kelvin to Celsius - # - # Check and correct for negative values - if (Kelvin < 0) { - # If Kelvin is less than zero, set it to zero - Kelvin <- 0 - } - Celsius <- Kelvin - 273.15 - return(Celsius) -} - - -################################################## -################# as.radians ##################### -################################################## -as.radians <- function(degrees) { - # Converts from degrees to radians - radians <- degrees * (pi / 180) - return(radians) -} - - - -################################################## -################# as.degrees ##################### -################################################## -as.degrees <- function(radians) { - # Converts from radians to degrees - degrees <- radians * (180 / pi) - return(degrees) -} - - - -################################################## -############### molarity2mass #################### -################################################## -molarity2mass <- function(formulamass, volume, molarity) { - # Calculates the required mass of - # the substance to be dissolved. - # ARGS: formulamass - formula mass of the substance (in gram per mole) - # volume - volume of the final solution (in liters) - # molarity - molarity (in moles per liter) - # VALUE: mass of substance (in grams) - # - mass <- formulamass * volume * molarity - # Unit check: - # [g * mol-1] * [liter] * [mole * liter-1] = [g] - return(mass) -} - - - -################################################## -#################### AVS2SHE ##################### -################################################## -AVS2SHE <- function(avs) { - # Converts from absolute vacuum scale (AVS) to SHE scale - she <- -(4.5 + avs) - return(she) -} - - - -################################################## -#################### SHE2AVS ##################### -################################################## -SHE2AVS <- function(she) { - # Converts from SHE scale to absolute vacuum (AVS) scale - avs <- -(4.5 + she) - return(avs) -} - - - -################################################## -############### ConvertRefPotEC ################## -################################################## -ConvertRefPotEC <- function(argpotential, argrefscale, valuerefscale) { - # Converts from an electrochemical reference potential scale into another - # SHE: standard hydrogen electrode scale - # Ag/AgCl: silver silver-chloride electrode scale - # SCE: standard calomel scale - # - - ##### Add more reference electrodes here >> - refpotatSHEzero <- c( 0, -0.21, -0.24, 3) - refrownames <- c( "SHE", "Ag/AgCl", "SCE", "Li/Li+") - refcolnames <- c("SHE0", "AgCl0", "SCE0", "Li0") - ##### Add more reference electrodes here << - # - SHE0 <- data.frame(matrix(refpotatSHEzero, ncol=length(refpotatSHEzero), byrow=T)) - refpotmtx <- matrix(NA, length(SHE0), length(SHE0)) - refpotmtx[,1] <- matrix(as.matrix(SHE0), ncol=1, byrow=T) - for (c in 2:length(SHE0)) { - # loop over columns (except the first) - for (r in 1:length(SHE0)) { - # loop over rows - refpotmtx[r, c] <- refpotmtx[r, 1] - refpotmtx[c, 1] - } - } - refpotdf <- as.data.frame(refpotmtx) - names(refpotdf) <- refcolnames - row.names(refpotdf) <- refrownames - ## So far we have made a matrix of all the possible combinations, - ## given the vector refpotatSHEzero. The matrix is not strictly necessary, - ## but it may prove useful later. It does. - # - # Match argrefscale to the refrownames - argmatch <- match(argrefscale, refrownames, nomatch = 0) - # Match valuerefscale to the refrownames - valuematch <- match(valuerefscale, refrownames, nomatch = 0) - # We simply assume that the match was well-behaved - valuepotential <- argpotential + refpotdf[valuematch, argmatch] - # Check that arg and value electrodes are within bounds for a match - if (argmatch == 0 || valuematch == 0) { - # No match - # Perform suitable action - message("Arg out of bounds in call to ConvertRefPot") - valuepotential <- NA - } - return(valuepotential) -} - - -################################################## -################# ConvertRefPot ################## -################################################## -ConvertRefPot <- function(argpotential, argrefscale, valuerefscale) { - # Check that argpotential is valid numeric - - # IDEA: make a matrix out of these (scale names and flags) - - # Valid scales - scale.names <- list() - scale.names[["SHE"]] <- c("SHE", "NHE", "she", "nhe") - scale.names[["AgCl"]] <- c("Ag/AgCl", "AgCl", "ag/agcl", "agcl") - scale.names[["SCE"]] <- c("SCE", "sce") - scale.names[["Li"]] <- c("Li/Li+", "Li", "Li+", "li", "li+", "li/li+") - scale.names[["AVS"]] <- c("AVS", "avs") - - # Set flags - bool.flags <- as.data.frame(matrix(0, nrow = length(scale.names), ncol = 2)) - names(bool.flags) <- c("argref", "valueref") - row.names(bool.flags) <- names(scale.names) - - # argrefscale - # Check that argrefscale is valid character mode - # ... - - for (j in 1:length(row.names(bool.flags))) { - if (any(scale.names[[row.names(bool.flags)[j]]] == argrefscale)) { - bool.flags[row.names(bool.flags)[j], "argref"] <- j - } - } - - - # valuerefscale - # Check that valuerefscale is valid character mode - # ... - - for (k in 1:length(row.names(bool.flags))) { - if (any(scale.names[[row.names(bool.flags)[k]]] == valuerefscale)) { - bool.flags[row.names(bool.flags)[k], "valueref"] <- k - } - } - - # Depending on which flags are set, call the corresponding function - - decision.vector <- colSums(bool.flags) - - # Check if both scales are the same (no conversion needed). If so, abort gracefully. - # ... - - if (decision.vector["argref"] == 5 || decision.vector["valueref"] == 5) { - # AVS is requested, deal with it it - if (decision.vector["argref"] == 5) { - # Conversion _from_ AVS - rnpotential <- ConvertRefPotEC(AVS2SHE(argpotential), - "SHE", - scale.names[[decision.vector["valueref"]]][1]) - } - if (decision.vector["valueref"] == 5) { - # Conversion _to_ AVS - rnpotential <- SHE2AVS(ConvertRefPotEC(argpotential, - scale.names[[decision.vector["argref"]]][1], - "SHE")) - } - } else { - rnpotential <- ConvertRefPotEC(argpotential, - scale.names[[decision.vector["argref"]]][1], - scale.names[[decision.vector["valueref"]]][1]) - } - return(rnpotential) -} \ No newline at end of file diff --git a/common/AVS2SHE.R b/common/AVS2SHE.R new file mode 100644 index 0000000..6fc2d06 --- /dev/null +++ b/common/AVS2SHE.R @@ -0,0 +1,8 @@ +################################################## +#################### AVS2SHE ##################### +################################################## +AVS2SHE <- function(avs) { + # Converts from absolute vacuum scale (AVS) to SHE scale + she <- -(4.5 + avs) + return(she) +} \ No newline at end of file diff --git a/common/Celsius2Kelvin.R b/common/Celsius2Kelvin.R new file mode 100644 index 0000000..4c06661 --- /dev/null +++ b/common/Celsius2Kelvin.R @@ -0,0 +1,14 @@ +################################################## +############### Celsius2Kelvin ################### +################################################## +Celsius2Kelvin <- function(Celsius) { + # Converts temperature from Celsius to Kelvin + # + # Check and correct for values below -273.15 + if (Celsius < -273.15) { + # If Celsis is less than absolute zero, set it to absolute zero + Celsius <- -273.15 + } + Kelvin <- Celsius + 273.15 + return(Kelvin) +} \ No newline at end of file diff --git a/common/ConvertRefPot.R b/common/ConvertRefPot.R new file mode 100644 index 0000000..1143d63 --- /dev/null +++ b/common/ConvertRefPot.R @@ -0,0 +1,74 @@ +source(SHE2AVS.R) +source(AVS2SHE.R) +source(ConvertRefPotEC.R) + +################################################## +################# ConvertRefPot ################## +################################################## +ConvertRefPot <- function(argpotential, argrefscale, valuerefscale) { + # Check that argpotential is valid numeric + + # IDEA: make a matrix out of these (scale names and flags) + + # Valid scales + scale.names <- list() + scale.names[["SHE"]] <- c("SHE", "NHE", "she", "nhe") + scale.names[["AgCl"]] <- c("Ag/AgCl", "AgCl", "ag/agcl", "agcl") + scale.names[["SCE"]] <- c("SCE", "sce") + scale.names[["Li"]] <- c("Li/Li+", "Li", "Li+", "li", "li+", "li/li+") + scale.names[["AVS"]] <- c("AVS", "avs") + + # Set flags + bool.flags <- as.data.frame(matrix(0, nrow = length(scale.names), ncol = 2)) + names(bool.flags) <- c("argref", "valueref") + row.names(bool.flags) <- names(scale.names) + + # argrefscale + # Check that argrefscale is valid character mode + # ... + + for (j in 1:length(row.names(bool.flags))) { + if (any(scale.names[[row.names(bool.flags)[j]]] == argrefscale)) { + bool.flags[row.names(bool.flags)[j], "argref"] <- j + } + } + + + # valuerefscale + # Check that valuerefscale is valid character mode + # ... + + for (k in 1:length(row.names(bool.flags))) { + if (any(scale.names[[row.names(bool.flags)[k]]] == valuerefscale)) { + bool.flags[row.names(bool.flags)[k], "valueref"] <- k + } + } + + # Depending on which flags are set, call the corresponding function + + decision.vector <- colSums(bool.flags) + + # Check if both scales are the same (no conversion needed). If so, abort gracefully. + # ... + + if (decision.vector["argref"] == 5 || decision.vector["valueref"] == 5) { + # AVS is requested, deal with it it + if (decision.vector["argref"] == 5) { + # Conversion _from_ AVS + rnpotential <- ConvertRefPotEC(AVS2SHE(argpotential), + "SHE", + scale.names[[decision.vector["valueref"]]][1]) + } + if (decision.vector["valueref"] == 5) { + # Conversion _to_ AVS + rnpotential <- SHE2AVS(ConvertRefPotEC(argpotential, + scale.names[[decision.vector["argref"]]][1], + "SHE")) + } + } else { + rnpotential <- ConvertRefPotEC(argpotential, + scale.names[[decision.vector["argref"]]][1], + scale.names[[decision.vector["valueref"]]][1]) + } + return(rnpotential) +} \ No newline at end of file diff --git a/common/ConvertRefPotEC.R b/common/ConvertRefPotEC.R new file mode 100644 index 0000000..4b4ec54 --- /dev/null +++ b/common/ConvertRefPotEC.R @@ -0,0 +1,48 @@ +################################################## +############### ConvertRefPotEC ################## +################################################## +ConvertRefPotEC <- function(argpotential, argrefscale, valuerefscale) { + # Converts from an electrochemical reference potential scale into another + # SHE: standard hydrogen electrode scale + # Ag/AgCl: silver silver-chloride electrode scale + # SCE: standard calomel scale + # + + ##### Add more reference electrodes here >> + refpotatSHEzero <- c( 0, -0.21, -0.24, 3) + refrownames <- c( "SHE", "Ag/AgCl", "SCE", "Li/Li+") + refcolnames <- c("SHE0", "AgCl0", "SCE0", "Li0") + ##### Add more reference electrodes here << + # + SHE0 <- data.frame(matrix(refpotatSHEzero, ncol=length(refpotatSHEzero), byrow=T)) + refpotmtx <- matrix(NA, length(SHE0), length(SHE0)) + refpotmtx[,1] <- matrix(as.matrix(SHE0), ncol=1, byrow=T) + for (c in 2:length(SHE0)) { + # loop over columns (except the first) + for (r in 1:length(SHE0)) { + # loop over rows + refpotmtx[r, c] <- refpotmtx[r, 1] - refpotmtx[c, 1] + } + } + refpotdf <- as.data.frame(refpotmtx) + names(refpotdf) <- refcolnames + row.names(refpotdf) <- refrownames + ## So far we have made a matrix of all the possible combinations, + ## given the vector refpotatSHEzero. The matrix is not strictly necessary, + ## but it may prove useful later. It does. + # + # Match argrefscale to the refrownames + argmatch <- match(argrefscale, refrownames, nomatch = 0) + # Match valuerefscale to the refrownames + valuematch <- match(valuerefscale, refrownames, nomatch = 0) + # We simply assume that the match was well-behaved + valuepotential <- argpotential + refpotdf[valuematch, argmatch] + # Check that arg and value electrodes are within bounds for a match + if (argmatch == 0 || valuematch == 0) { + # No match + # Perform suitable action + message("Arg out of bounds in call to ConvertRefPot") + valuepotential <- NA + } + return(valuepotential) +} \ No newline at end of file diff --git a/common/Kelvin2Celsius.R b/common/Kelvin2Celsius.R new file mode 100644 index 0000000..2c688b4 --- /dev/null +++ b/common/Kelvin2Celsius.R @@ -0,0 +1,14 @@ +################################################## +############### Kelvin2Celsius ################### +################################################## +Kelvin2Celsius <- function(Kelvin) { + # Converts temperature from Kelvin to Celsius + # + # Check and correct for negative values + if (Kelvin < 0) { + # If Kelvin is less than zero, set it to zero + Kelvin <- 0 + } + Celsius <- Kelvin - 273.15 + return(Celsius) +} \ No newline at end of file diff --git a/common/ProvideSampleId.R b/common/ProvideSampleId.R new file mode 100644 index 0000000..6b8b7c1 --- /dev/null +++ b/common/ProvideSampleId.R @@ -0,0 +1,25 @@ +################################################## +################ ProvideSampleId ################# +################################################## +ProvideSampleId <- function (fullpathwithfilename) { + ### OBS! Only very rudimentary error-checking. + ### If the filename is formatted as \w*-\w*-\w*, we use the middle segment, + ### otherwise we use the whole string (excluding the extension) + # Extract the name of the parent directory of the datafilename argument + substrateid <- basename(dirname(fullpathwithfilename)) + # Extract the name of the method from the filename-part + # First split the filename over all hyphens + nameparts <- strsplit(basename(fullpathwithfilename), "-")[[1]] + # If the number of nameparts exceed 3, save the whole filename as methodid, otherwise use the middle part + if (length(nameparts) > 3) { + # We need to lose the file extension from the last namepart + nameparts[length(nameparts)] <- strsplit(nameparts[length(nameparts)], "\\.")[[1]][1] + methodid <- paste(nameparts, collapse = "-") + } else { + methodid <- nameparts[2] + } + # Make an informative sampleid + sampleid <- paste(substrateid, methodid, sep = "-") + # + return(sampleid) +} \ No newline at end of file diff --git a/common/SHE2AVS.R b/common/SHE2AVS.R new file mode 100644 index 0000000..31c7e48 --- /dev/null +++ b/common/SHE2AVS.R @@ -0,0 +1,8 @@ +################################################## +#################### SHE2AVS ##################### +################################################## +SHE2AVS <- function(she) { + # Converts from SHE scale to absolute vacuum (AVS) scale + avs <- -(4.5 + she) + return(avs) +} \ No newline at end of file diff --git a/common/as.degrees.R b/common/as.degrees.R new file mode 100644 index 0000000..3a80188 --- /dev/null +++ b/common/as.degrees.R @@ -0,0 +1,8 @@ +################################################## +################# as.degrees ##################### +################################################## +as.degrees <- function(radians) { + # Converts from radians to degrees + degrees <- radians * (180 / pi) + return(degrees) +} \ No newline at end of file diff --git a/common/as.radians.R b/common/as.radians.R new file mode 100644 index 0000000..7f643cc --- /dev/null +++ b/common/as.radians.R @@ -0,0 +1,8 @@ +################################################## +################# as.radians ##################### +################################################## +as.radians <- function(degrees) { + # Converts from degrees to radians + radians <- degrees * (pi / 180) + return(radians) +} \ No newline at end of file diff --git a/common/common.info b/common/common.info new file mode 100644 index 0000000..2f8cf0b --- /dev/null +++ b/common/common.info @@ -0,0 +1,19 @@ +# 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 diff --git a/common/int2padstr.R b/common/int2padstr.R new file mode 100644 index 0000000..aa9aba4 --- /dev/null +++ b/common/int2padstr.R @@ -0,0 +1,21 @@ +################################################## +################## int2padstr #################### +################################################## +int2padstr <- function (ii, pchr, w) { + ## Description: + ## Converts an integer or a vector of integers to + ## a string padded with characters. + ## Usage: + ## int2padstr(ii, pchr, w) + ## Arguments: + ## ii: integer or vector of integers + ## pchr: a padding character (e.g., "0") + ## w: width of the return string (an integer) + ## Make sure to set the width longer than + ## or equal to the length of the biggest integer. + ## For example, if the integers (ii) are + ## in the range 1 - 100, set w to at least 3. + ## Value: + ## A character string or a vector of character strings + gsub(" ", pchr, formatC(ii, format="s", mode="character", width = w)) +} \ No newline at end of file diff --git a/common/molarity2mass.R b/common/molarity2mass.R new file mode 100644 index 0000000..5eb8d49 --- /dev/null +++ b/common/molarity2mass.R @@ -0,0 +1,16 @@ +################################################## +############### molarity2mass #################### +################################################## +molarity2mass <- function(formulamass, volume, molarity) { + # Calculates the required mass of + # the substance to be dissolved. + # ARGS: formulamass - formula mass of the substance (in gram per mole) + # volume - volume of the final solution (in liters) + # molarity - molarity (in moles per liter) + # VALUE: mass of substance (in grams) + # + mass <- formulamass * volume * molarity + # Unit check: + # [g * mol-1] * [liter] * [mole * liter-1] = [g] + return(mass) +} \ No newline at end of file diff --git a/common/unborn/It2charge.R b/common/unborn/It2charge.R new file mode 100644 index 0000000..6ef135b --- /dev/null +++ b/common/unborn/It2charge.R @@ -0,0 +1,35 @@ +################################################## +################### 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) +} \ No newline at end of file diff --git a/common/unborn/LinearBaseline.R b/common/unborn/LinearBaseline.R new file mode 100644 index 0000000..5b7f5a9 --- /dev/null +++ b/common/unborn/LinearBaseline.R @@ -0,0 +1,23 @@ +################################################## +################ 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) +} \ No newline at end of file diff --git a/xrdtf.R b/xrdtf.R deleted file mode 100644 index bb78877..0000000 --- a/xrdtf.R +++ /dev/null @@ -1,808 +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 - - - - -################################################## -################## 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) -} - - - - -################################################## -################## 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(), source("/home/taha/chepec/chetex/common/R/common.R") - D <- (shapeconstant * wavelength) / (as.radians(integralbreadth) * cos(as.radians(thth))) - # cos() - angles must be in radians, not degrees! - return(D) #units of angstrom -} - - - - -################################################## -################### pdf2df ####################### -################################################## -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, this function simply strips the less-than character. - # (Perhaps not true, see the int.Tex column) - # 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) - # attr: This function sets the following attributes: - # ApplicationName, - # ApplicationVersion, - # pdfNumber, - # 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 = "") - )) - } - # - 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) -} - - - - -################################################## -################### 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 -} - - - - - -################################################## -#################### uxd2df ###################### -################################################## -uxd2df <- function(uxdfile) { - # Function for reading UXD files # Assumptions: data in two columns - # Args: uxdfile (filename with extension) - # Returns: dataframe with three 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] - - rgxp.sampleid <- "[^/]*(?=\\.\\w*)" ## THIS REQUIRES perl=TRUE - # Regular expression that extracts the filename out of a full path. - # Matches and extracts everything from the last forward slash (assuming Unix slashes) - # up until a dot folllowed by an arbitrary number of alphanumeric characters. - sampleidmtch <- regexpr(rgxp.sampleid, uxdfile, perl=TRUE) - # Check that there was a match - if (sampleidmtch < 0) { - # -1 means no match - sampleid <- uxdfile - # If match was unsuccessful we use the argument as passed to this function as sampleid - } - sampleid <- substr(uxdfile, sampleidmtch, (sampleidmtch + attr(sampleidmtch, "match.length") - 1)) - - zz <- textConnection(f, "r") - ff <- data.frame(sampleid, matrix(scan(zz, - what = numeric()), ncol=2, byrow=T)) - names(ff) <- c("sampleid", "angle", "intensity") - close(zz) - - #zz <- file(datafile, "w") #open connection to datafile - #write.table(ff, file=datafile, row.names=F, sep=",") - #close(zz) - - # Return dataframe - ff -} - - - - -################################################## -################### muxd2df ###################### -################################################## -muxd2df <- function(uxdfile, range.descriptor) { - # Function that reads an UXD file which contains several ranges - # (created in a programmed run, for example) - # Arguments - # :: uxdfile (filename with extension) - # :: range.descriptor (an array with as many elements as - # there are ranges in the uxdfile) - # Returns: dataframe with 3 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) - - # 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)) { - zz <- textConnection(f[starts[s]:ends[s]], "r") - ff <- rbind(ff, data.frame(range.descriptor[s], - matrix(scan(zz, what = numeric()), ncol=2, byrow=T))) - close(zz) - } - names(ff) <- c("sampleid", "angle", "intensity") - - # Return dataframe - ff -} - - - - -################################################## -################### 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 -} - - - -################################################## -################### 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) -} - - - - - - - - -# -------- ##################### -------- # -# -------- #### REPAIR SHOP #### -------- # -# -------- ##################### -------- # - - - -################################################## -################ 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 #### -} - - - - -################################################## -################ 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 -} - - - -################################################## -################## 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) -} - - - - -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 "))) -} - - -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 -}