R-common/XRD-TF/unborn/EliminateKa2()

##################################################
################ 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 ####
}