Reading the absorbance data from text file, and linear/spline interpolation function:
absorbance <- read.table("absorbance.txt", head = TRUE)
absorbance <- rbind(c(200, 0), c(300, 0), absorbance, c(1000, 0), c(1100, 0)) # make sure it goes to 0 on each side
# absorption <- approxfun(1e7/absorbance$x, absorbance$y, yleft=0,
# yright=0)
absorption <- splinefun(1e+07/absorbance$x, absorbance$y, method = "monoH.FC")
x <- seq(1e+07/700, 1e+07/500, length = 200)
x2 <- 1e+07/x
y <- absorption(x)
plot(x, y, t = "l", xlab = expression("Wavenumber / " * cm^-1), ylab = "Absorbance")
Now we define the function phi, and use it to get the REP for both modes
phi <- function(nu, epsilon = 5, maxi = 21000) {
fun <- function(x) absorption(x)/(x * (x - nu))
I1 <- integrate(fun, epsilon, nu - epsilon, rel.tol = 1e-08, subdivisions = 200L)
I2 <- integrate(fun, nu + epsilon, maxi, rel.tol = 1e-08, subdivisions = 200L)
I1$value + I2$value - (0+1i) * pi * absorption(nu)/nu
}
rep <- function(deltanu = 595) {
phi1 <- sapply(x, phi)
phi2 <- sapply(x - deltanu, phi)
(x - deltanu)^4 * Mod(phi1 - phi2)^2
}
rep1 <- rep(595)
rep2 <- rep(1650)
plot(x2, rep1/max(rep1), t = "l", col = 2, xlab = "Wavelength / nm", ylab = "Normalised intensity")
lines(x2, rep2/max(rep2), col = 3)
lines(x2, absorption(x)/max(absorption(x)), lty = 2)
legend("topleft", lty = c(2, 1, 1), col = c(1, 2, 3), c("Absorbance", "595 mode",
"1650 mode"))
