# d <- read.table('ph5004237_si_002.txt')
d2 <- read.table('Olmon_PRB2012_SC.dat', skip=3)
wat <- read.table('water.txt')
gold <- data.frame(wavelength = d$V1, imag = d$V5, real = d$V4)
m <- subset(melt(gold, id="wavelength"), wavelength < 1000)
gold2 <- data.frame(wavelength = d2$V2*1e9, imag = d2$V4, real = d2$V3)
m2 <- subset(melt(gold2, id="wavelength"), wavelength > min(m$wavelength) &
wavelength < max(m$wavelength))
library(ggplot2)
library(reshape2)
ggplot(m, aes(wavelength, value, colour=variable)) +
geom_line() +
geom_line(data=m2, linetype=2)
epsWater <- function(wavelength){
raw <- structure(list(wavelength = 300:800, V1 = c(1.8474, 1.8468, 1.8461,
1.8454, 1.8447, 1.8441, 1.8434, 1.8428, 1.8422, 1.8415, 1.8409,
1.8403, 1.8397, 1.8391, 1.8385, 1.8379, 1.8373, 1.8368, 1.8362,
1.8357, 1.8351, 1.8346, 1.834, 1.8335, 1.833, 1.8324, 1.8319,
1.8314, 1.8309, 1.8304, 1.8299, 1.8294, 1.829, 1.8285, 1.828,
1.8275, 1.8271, 1.8266, 1.8262, 1.8257, 1.8253, 1.8248, 1.8244,
1.824, 1.8236, 1.8231, 1.8227, 1.8223, 1.8219, 1.8215, 1.8211,
1.8207, 1.8203, 1.8199, 1.8195, 1.8192, 1.8188, 1.8184, 1.818,
1.8177, 1.8173, 1.8169, 1.8166, 1.8162, 1.8159, 1.8155, 1.8152,
1.8148, 1.8145, 1.8142, 1.8138, 1.8135, 1.8132, 1.8129, 1.8125,
1.8122, 1.8119, 1.8116, 1.8113, 1.811, 1.8107, 1.8104, 1.8101,
1.8098, 1.8095, 1.8092, 1.8089, 1.8086, 1.8084, 1.8081, 1.8078,
1.8075, 1.8073, 1.807, 1.8067, 1.8064, 1.8062, 1.8059, 1.8057,
1.8054, 1.8051, 1.8049, 1.8046, 1.8044, 1.8041, 1.8039, 1.8037,
1.8034, 1.8032, 1.8029, 1.8027, 1.8025, 1.8022, 1.802, 1.8018,
1.8015, 1.8013, 1.8011, 1.8009, 1.8006, 1.8004, 1.8002, 1.8,
1.7998, 1.7996, 1.7994, 1.7991, 1.7989, 1.7987, 1.7985, 1.7983,
1.7981, 1.7979, 1.7977, 1.7975, 1.7973, 1.7971, 1.7969, 1.7968,
1.7966, 1.7964, 1.7962, 1.796, 1.7958, 1.7956, 1.7954, 1.7953,
1.7951, 1.7949, 1.7947, 1.7946, 1.7944, 1.7942, 1.794, 1.7939,
1.7937, 1.7935, 1.7933, 1.7932, 1.793, 1.7928, 1.7927, 1.7925,
1.7924, 1.7922, 1.792, 1.7919, 1.7917, 1.7916, 1.7914, 1.7913,
1.7911, 1.7909, 1.7908, 1.7906, 1.7905, 1.7903, 1.7902, 1.79,
1.7899, 1.7898, 1.7896, 1.7895, 1.7893, 1.7892, 1.789, 1.7889,
1.7888, 1.7886, 1.7885, 1.7884, 1.7882, 1.7881, 1.7879, 1.7878,
1.7877, 1.7875, 1.7874, 1.7873, 1.7872, 1.787, 1.7869, 1.7868,
1.7866, 1.7865, 1.7864, 1.7863, 1.7861, 1.786, 1.7859, 1.7858,
1.7856, 1.7855, 1.7854, 1.7853, 1.7852, 1.7851, 1.7849, 1.7848,
1.7847, 1.7846, 1.7845, 1.7844, 1.7842, 1.7841, 1.784, 1.7839,
1.7838, 1.7837, 1.7836, 1.7835, 1.7834, 1.7832, 1.7831, 1.783,
1.7829, 1.7828, 1.7827, 1.7826, 1.7825, 1.7824, 1.7823, 1.7822,
1.7821, 1.782, 1.7819, 1.7818, 1.7817, 1.7816, 1.7815, 1.7814,
1.7813, 1.7812, 1.7811, 1.781, 1.7809, 1.7808, 1.7807, 1.7806,
1.7805, 1.7804, 1.7803, 1.7802, 1.7801, 1.78, 1.7799, 1.7799,
1.7798, 1.7797, 1.7796, 1.7795, 1.7794, 1.7793, 1.7792, 1.7791,
1.779, 1.779, 1.7789, 1.7788, 1.7787, 1.7786, 1.7785, 1.7784,
1.7783, 1.7783, 1.7782, 1.7781, 1.778, 1.7779, 1.7778, 1.7778,
1.7777, 1.7776, 1.7775, 1.7774, 1.7774, 1.7773, 1.7772, 1.7771,
1.777, 1.777, 1.7769, 1.7768, 1.7767, 1.7766, 1.7766, 1.7765,
1.7764, 1.7763, 1.7762, 1.7762, 1.7761, 1.776, 1.7759, 1.7759,
1.7758, 1.7757, 1.7756, 1.7756, 1.7755, 1.7754, 1.7754, 1.7753,
1.7752, 1.7751, 1.7751, 1.775, 1.7749, 1.7748, 1.7748, 1.7747,
1.7746, 1.7746, 1.7745, 1.7744, 1.7744, 1.7743, 1.7742, 1.7741,
1.7741, 1.774, 1.7739, 1.7739, 1.7738, 1.7737, 1.7737, 1.7736,
1.7735, 1.7735, 1.7734, 1.7733, 1.7733, 1.7732, 1.7731, 1.7731,
1.773, 1.7729, 1.7729, 1.7728, 1.7728, 1.7727, 1.7726, 1.7726,
1.7725, 1.7724, 1.7724, 1.7723, 1.7722, 1.7722, 1.7721, 1.7721,
1.772, 1.7719, 1.7719, 1.7718, 1.7717, 1.7717, 1.7716, 1.7716,
1.7715, 1.7714, 1.7714, 1.7713, 1.7713, 1.7712, 1.7711, 1.7711,
1.771, 1.771, 1.7709, 1.7709, 1.7708, 1.7707, 1.7707, 1.7706,
1.7706, 1.7705, 1.7704, 1.7704, 1.7703, 1.7703, 1.7702, 1.7702,
1.7701, 1.7701, 1.77, 1.7699, 1.7699, 1.7698, 1.7698, 1.7697,
1.7697, 1.7696, 1.7696, 1.7695, 1.7694, 1.7694, 1.7693, 1.7693,
1.7692, 1.7692, 1.7691, 1.7691, 1.769, 1.769, 1.7689, 1.7688,
1.7688, 1.7687, 1.7687, 1.7686, 1.7686, 1.7685, 1.7685, 1.7684,
1.7684, 1.7683, 1.7683, 1.7682, 1.7682, 1.7681, 1.7681, 1.768,
1.768, 1.7679, 1.7679, 1.7678, 1.7678, 1.7677, 1.7677, 1.7676,
1.7676, 1.7675, 1.7675, 1.7674, 1.7674, 1.7673, 1.7673, 1.7672,
1.7672, 1.7671, 1.7671, 1.767, 1.767, 1.7669, 1.7669, 1.7668,
1.7668, 1.7667, 1.7667, 1.7666, 1.7666, 1.7665, 1.7665, 1.7664,
1.7664, 1.7663, 1.7663, 1.7662, 1.7662, 1.7661, 1.7661, 1.766,
1.766, 1.7659, 1.7659, 1.7659, 1.7658, 1.7658, 1.7657, 1.7657,
1.7656, 1.7656, 1.7655, 1.7655, 1.7654, 1.7654, 1.7653, 1.7653,
1.7652, 1.7652)), .Names = c("wavelength", "V1"), row.names = c(NA,
-501L), class = "data.frame")
data.frame(wavelength=wavelength, epsilon = approx(raw[,1],raw[,2],wavelength)$y)
}
epsGold <- function(wavelength){
epsr <- approx(gold$wavelength,gold$real,wavelength)$y
epsi <- approx(gold$wavelength,gold$imag,wavelength)$y
data.frame(wavelength=wavelength,
epsilon = epsr + 1i*epsi)
}
epsGold2 <- function(wavelength){
epsr <- approx(gold2$wavelength,gold2$real,wavelength)$y
epsi <- approx(gold2$wavelength,gold2$imag,wavelength)$y
data.frame(wavelength=wavelength,
epsilon = epsr + 1i*epsi)
}
library(mie)
library(dielectric)
goldJC <- epsAu(seq(400, 800))
goldACS <- epsGold(seq(400, 800))
goldPRB <- epsGold2(seq(400, 800))
water <- sqrt(epsWater(seq(400, 800)))[["epsilon"]]
xsec1 <- with(goldJC, mie(wavelength, epsilon, radius=29, medium=water, efficiency=FALSE))
xsec2 <- with(goldACS, mie(wavelength, epsilon, radius=29, medium=water, efficiency=FALSE))
xsec3 <- with(goldPRB, mie(wavelength, epsilon, radius=29, medium=water, efficiency=FALSE))
# matplot(cross_sections$wavelength, cross_sections[, -1], type="l", lty=1,
# xlab=expression(lambda/mu*m), ylab=expression(sigma/mu*m^2))
# matlines(cross_sections2$wavelength, cross_sections2[, -1], lty=2)
# matlines(cross_sections3$wavelength, cross_sections3[, -1], lty=3)
#
# legend("topright", names(cross_sections)[-1], col=1:3, lty=1)
# xsec1$name <- "JC"
xsec2$name <- "ACS"
xsec3$name <- "PRB"
all <- rbind(xsec2, xsec3)
mm <- melt(all, id=c("wavelength","name"))
ggplot(mm, aes(wavelength, value, group=interaction(variable,name))) +
geom_line(aes(colour=name)) +
facet_wrap(~variable, scales="free") +
theme_minimal() +
labs(y="",colour="")
ggsave("epsAu.pdf",width=12,height=3)
library(cda)
library(cubature)
library(egg)
library(tidyr)
library(purrr)
library(dplyr)
library(tibble)
library(dielectric)
weight <- function(a, sd=10){
dnorm(a, mean = 28, sd = sd)
}
tibble(a = seq(1,100,length=100)) %>%
mutate(p = weight(a)) %>%
ggplot(aes(a, p)) + geom_line()
## check that the area is 1
integrate(weight, -Inf, Inf)
#
# mie(goldACS$wavelength, goldACS$epsilon, radius=30,
# medium=water, efficiency=FALSE)
wrap_model <- function(p, material=goldACS, sd=5,...){
xsec1 <- mie(material$wavelength, material$epsilon, radius=p,
medium=water, efficiency=FALSE)
w <- weight(p,sd=sd)
c(xsec1$extinction, xsec1$absorption, xsec1$scattering)* w
}
length(wrap_model(30))
wrap_me <- function(sd){
tmp <- pcubature(wrap_model, sd=sd, material=goldPRB, lowerLimit = 1, upperLimit = 100,
tol=1e-3, maxEval = 300, norm = "PAIRED",
fDim = 3*nrow(goldACS) )
tibble(wavelength = rep(goldACS$wavelength, 3),
value = tmp$integral,
variable = rep(c("extinction", "absorption", "scattering"),
each=nrow(goldACS)),
name='distribution'
)
}
all <- plyr::mdply(data.frame(sd=seq(2,8,by=2)), wrap_me)
# quadrature <- wrap_me(5)
raw <- setNames(read.table("raw.txt"),
c("wavelength", "scattering", "absorption", "extinction"))
md <- melt(raw, id=1)
md$name <- "data"
ggplot(subset(mm, name=="PRB"), aes(wavelength, value, group=interaction(variable,name))) +
geom_line(data=md,aes(y=value*4e4),lwd=1.2,col="grey70") +
# geom_line(col='black') +
geom_line(data=all, aes(colour=factor(sd),group=sd),lty=2) +
facet_wrap(~variable, scales="free") +
theme_minimal() +
labs(y="",colour="sd")
# ggsave("epsAu-disp2.pdf",width=12,height=3)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.