Resonance Energy Profile

absorption <- function(x) {
    dnorm(x, mean = 1e+07/700, sd = 500)
}

x <- seq(1e+07/800, 1e+07/500, length = 200)
x2 <- 1e+07/x
y <- absorption(x)
abs <- data.frame(wavelength = x2, absorption = y)
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 as a function of Raman shift.

phi <- function(nu, epsilon = 10, mini = 5000, maxi = 21000) {

    fun <- function(x) absorption(x)/(x * (x - nu))

    I1 <- integrate(fun, mini, nu - epsilon, rel.tol = 1e-10, subdivisions = 10000L)
    I2 <- integrate(fun, nu + epsilon, maxi, rel.tol = 1e-10, subdivisions = 10000L)

    (I1$value + I2$value) - (0+1i) * pi * absorption(nu)/nu

}

rep <- function(deltanu = 595, x) {
    phi1 <- sapply(x, phi)
    phi2 <- sapply(x - deltanu, phi)
    data.frame(wavenumber = x, wavelength = x2, rep = (x - deltanu)^4 * Mod(phi1 - 
        phi2)^2)
}

params <- data.frame(deltanu = c(100, 500, 1000, 2000, 3000, 4000))
results <- mdply(params, rep, x = x)
results <- ddply(results, "deltanu", mutate, normalised = rep/max(rep))

ggplot(results, aes(wavelength, rep)) + geom_line(aes(colour = deltanu, group = deltanu))

ggplot(results, aes(wavelength, normalised)) + geom_line(aes(colour = deltanu, 
    group = deltanu)) + geom_line(aes(wavelength, absorption/max(absorption)), 
    data = abs, linetype = "dashed")

x2 <- seq(800, 500, length = 200)
x <- 1e+07/x2
params <- data.frame(deltanu = seq(10, 5000, length = 200))
results <- mdply(params, rep, x = x)
results <- ddply(results, "deltanu", mutate, normalised = rep/max(rep))

ggplot(results, aes(wavelength, deltanu)) + geom_tile(aes(fill = normalised))