From 74450ff30c553c16d38bf1db03ff363b113d123b Mon Sep 17 00:00:00 2001 From: Taha Ahmed Date: Fri, 4 Mar 2011 11:07:31 +0100 Subject: [PATCH] Most recent change: added working electrode area to electrochemical functions. --- .gitignore | 3 + CHI.R | 20 +++- Renishaw.R | 2 + common.R | 6 +- xrdtf.R | 273 +++++++++++++++++++++++++++++++++++++++++++++-------- 5 files changed, 261 insertions(+), 43 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5d05fb6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.RData +*.Rhistory +*.Rhistory.save diff --git a/CHI.R b/CHI.R index d94b848..5e13ae6 100644 --- a/CHI.R +++ b/CHI.R @@ -158,7 +158,7 @@ chronocm2df <- function(datafilename) { ################################################## ################# chronoamp2df ################### ################################################## -chronoamp2df <- function(datafilename) { +chronoamp2df <- function(datafilename, wearea = 1) { # Function description: chronoamperometry data # CH Instruments potentiostat records all data using standard SI units, # so all potential values are in volts, currents are in amperes, @@ -201,6 +201,9 @@ chronoamp2df <- function(datafilename) { close(zz) } names(ff) <- c("step", "time", "current") + # Calculate current density + currentdensity <- ff$current / wearea + ff <- cbind(ff, currentdensity = currentdensity) # ### Collect attributes of this experiment # These attributes are specific for each kind of experiment, @@ -239,7 +242,7 @@ chronoamp2df <- function(datafilename) { ################################################## ############### amperometry2df ################### ################################################## -amperometry2df <- function(datafilename) { +amperometry2df <- function(datafilename, wearea = 1) { # Function description: for recorded amperometric i-T curves # CH Instruments potentiostat records all data using standard SI units, # so all potential values are in volts, currents are in amperes, @@ -281,6 +284,9 @@ amperometry2df <- function(datafilename) { close(zz) } names(ff) <- c("time", "current") + # Calculate current density + currentdensity <- ff$current / wearea + ff <- cbind(ff, currentdensity = currentdensity) # ### Collect attributes of this experiment # These attributes are specific for each kind of experiment, @@ -312,7 +318,7 @@ amperometry2df <- function(datafilename) { ################################################## #################### cv2df ####################### ################################################## -cv2df <- function(cvfilename) { +cv2df <- function(cvfilename, 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, @@ -355,6 +361,9 @@ cv2df <- function(cvfilename) { close(zz) } names(ff) <- c("segment", "cycle", "potential", "current", "charge") + # Calculate current density + currentdensity <- ff$current / wearea + ff <- cbind(ff[, 1:4], currentdensity = currentdensity, ff[, 5]) # ### Collect attributes of this experiment # These attributes are specific for each kind of experiment, @@ -394,7 +403,7 @@ cv2df <- function(cvfilename) { ################################################## ################### lsv2df ####################### ################################################## -lsv2df <- function(lsvfilename) { +lsv2df <- function(lsvfilename, 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, @@ -437,6 +446,9 @@ lsv2df <- function(lsvfilename) { close(zz) } names(ff) <- c("segment", "potential", "current", "charge") + # Calculate current density + currentdensity <- ff$current / wearea + ff <- cbind(ff[, 1:3], currentdensity = currentdensity, ff[, 4]) # ### Collect attributes of this experiment # These attributes are specific for each kind of experiment, diff --git a/Renishaw.R b/Renishaw.R index 5f9526b..f096ec3 100644 --- a/Renishaw.R +++ b/Renishaw.R @@ -22,6 +22,8 @@ Raman2df <- function(datafilename) { ##### # A nice algorithm that extracts the filename from the datafilename argument # and uses that as a sampleid in the returned dataframe + # THIS SHOULD PROBABLY BE CONVERTED INTO A STAND-ALONE FUNCTION + # Also make sure it works for vectors as well as single strings ##### rgxp.sampleid <- "[^/]*(?=\\.\\w*)" ## THIS REQUIRES perl=TRUE # Regular expression that extracts the filename out of a full path. diff --git a/common.R b/common.R index 0595261..7414748 100755 --- a/common.R +++ b/common.R @@ -74,6 +74,7 @@ Celsius2Kelvin <- function(Celsius) { Celsius <- -273.15 } Kelvin <- Celsius + 273.15 + return(Kelvin) } @@ -89,6 +90,7 @@ Kelvin2Celsius <- function(Kelvin) { Kelvin <- 0 } Celsius <- Kelvin - 273.15 + return(Celsius) } @@ -98,6 +100,7 @@ Kelvin2Celsius <- function(Kelvin) { as.radians <- function(degrees) { # Converts from degrees to radians radians <- degrees * (pi / 180) + return(radians) } @@ -106,7 +109,8 @@ as.radians <- function(degrees) { ################################################## as.degrees <- function(radians) { # Converts from radians to degrees - radians <- radians * (180 / pi) + degrees <- radians * (180 / pi) + return(degrees) } diff --git a/xrdtf.R b/xrdtf.R index b2e7e7a..6aced3c 100755 --- a/xrdtf.R +++ b/xrdtf.R @@ -11,6 +11,7 @@ # >>>> muxd2mtx # >>>> muxd2ls # - REPAIR SHOP +# - EliminateKa2 # - print.xtable.booktabs # - split.muxd # - strip.ka2 @@ -22,43 +23,6 @@ -################################################## -################ EliminateKa2 #################### -################################################## -EliminateKa2 <- function(xrdata) { - ##### STILL UNDER CONSTRUCTION #### - ##### STILL UNDER CONSTRUCTION #### - - # 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) - - ##### STILL UNDER CONSTRUCTION #### -} - - - - - ################################################## ################## matchpdf ###################### ################################################## @@ -321,8 +285,21 @@ uxd2df <- function(uxdfile) { 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(uxdfile, matrix(scan(zz, + ff <- data.frame(sampleid, matrix(scan(zz, what = numeric()), ncol=2, byrow=T)) names(ff) <- c("sampleid", "angle", "intensity") close(zz) @@ -507,6 +484,226 @@ muxd2ls <- function(uxdfile) { # -------- ##################### -------- # + +################################################## +################ 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 #################### ##################################################