Nothing
##' Multilayer Fresnel coefficients
##'
##' solves the EM problem of a multilayered interface
##' @title multilayer
##' @export
##' @param wavelength [vector] wavelength in nm
##' @param k0 [vector] wavevector in nm^-1
##' @param angle [vector] incident angles in radians
##' @param q [vector] normalised incident in-plane wavevector
##' @param epsilon list of N+2 dielectric functions, each of length 1 or length(wavelength)
##' @param thickness vector of N+2 layer thicknesses, first and last are dummy
##' @param d vector of distances where LFIEF are evaluated from each interface
##' @param dout vector of distances where LFIEF are evaluated outside the stack
##' @param polarisation [character] switch between p- and s- polarisation
##' @return fresnel coefficients and field profiles
##' @param ... unused
##' @author baptiste Auguie
##' @references
##' Principles of surface-enhanced Raman spectroscopy and related plasmonic effects.
##' Eric C. Le Ru and Pablo G. Etchegoin, published by Elsevier, Amsterdam (2009).
multilayer <- function(wavelength = 2*pi/k0, k0 = 2*pi/wavelength,
angle = asin(q), q = sin(angle),
epsilon = list(incident=1.5^2, 1.33),
thickness = c(0, 0),
polarisation = c('p', 's'),
d = 1, dout=d, ...){
## checks
stopifnot(thickness[1]==0L, thickness[length(thickness)]==0L)
polarisation <- match.arg(polarisation)
epsilon = do.call(cbind, epsilon)
## case pure scalars
if(nrow(epsilon) == 1L)
epsilon <- matrix(epsilon, nrow=length(k0),
ncol=length(thickness), byrow=TRUE)
## define constants
Nlambda <- length(k0)
Nq <- length(q)
Nlayer <- length(thickness)
k02 <- k0^2
kx <- outer(k0*sqrt(epsilon[,1] + 0i), q) # kx = q*k0
kx2 <- kx^2
## loop to calculate kiz
kiz <- array(0 + 0i, dim=c(Nlambda, Nq, Nlayer))
for (ii in seq(1, Nlayer)){
kiz[ , , ii] <- sqrt(matrix(epsilon[,ii]*k02, nrow=Nlambda, ncol=Nq) - kx2 + 0i)
}
## calculate the transition matrix M
M11 <- M22 <- 1 + 0i
M21 <- M12 <- 0 + 0i
Mi11 <- Mi12 <- Mi21 <- Mi22 <- array(1 + 0i, dim=c(Nlambda, Nq, Nlayer-1))
for (ii in seq(1, Nlayer-1)){
if(polarisation == 'p'){
Ki <- matrix(epsilon[,ii] / epsilon[,ii+1] + 0i , nrow=Nlambda, ncol=Nq) *
kiz[,,ii+1] / kiz[,,ii]
} else { # s-polarisation
Ki <- kiz[,,ii+1] / kiz[,,ii]
}
phasei <- exp(1i*thickness[ii]*kiz[,,ii])
Mi11[,,ii] <- 0.5*(1+Ki) / phasei
Mi21[,,ii] <- 0.5*(1-Ki) * phasei
Mi12[,,ii] <- 0.5*(1-Ki) / phasei
Mi22[,,ii] <- 0.5*(1+Ki) * phasei
M11new <- M11*Mi11[,,ii] + M12*Mi21[,,ii]
M21new <- M21*Mi11[,,ii] + M22*Mi21[,,ii]
M12new <- M11*Mi12[,,ii] + M12*Mi22[,,ii]
M22new <- M21*Mi12[,,ii] + M22*Mi22[,,ii]
M11 <- M11new
M12 <- M12new
M21 <- M21new
M22 <- M22new
}
## calculate the Fresnel coefficients
transmission <- 1 / M11
reflection <- M21 * transmission
## make a list of physical locations to sample the fields
Nd <- length(d)
sampling <- lapply(thickness, function(t) {
d[d<=t]
})
sampling[[1]] <- sampling[[length(sampling)]] <- dout # outside of the stack
## enhancement factors, default to 0
sizes <- lapply(sampling, length)
Ml.perp <- Ml.par <-
lapply(sizes[-length(sizes)], function(s) array(0, dim=c(Nlambda, Nq, s)))
Mr.perp <- Mr.par <-
lapply(sizes[-1], function(s) array(0, dim=c(Nlambda, Nq, s)))
## absolute positions, first interface is at 0
interfaces <- cumsum(thickness)[-length(thickness)]
distance <- list()
distance[[1]] <- -dout
#####################
## TM polarisation ##
#####################
if(polarisation == 'p'){
## field amplitudes p is for prime, reference is field in region 1
Hiy.H1y <- Hpiy.H1y <- Eix.E1 <- Epix.E1 <- Eiz.E1 <- Epiz.E1 <-
array(0+0i, dim=c(Nlambda, Nq, Nlayer))
Hiy.H1y[,,Nlayer] <- transmission
Hpiy.H1y[,,Nlayer] <- 0i
AuxE1 <- matrix(sqrt(epsilon[,1] + 0i) / k0 / epsilon[,Nlayer], nrow=Nlambda, ncol=Nq)
Eix.E1[,,Nlayer] <- Hiy.H1y[,,Nlayer] * kiz[,,Nlayer] * AuxE1
Epix.E1[,,Nlayer] <- 0i
AuxE2 <- outer(epsilon[,1] / epsilon[,Nlayer], Re(q))
Eiz.E1[,,Nlayer] <- - Hiy.H1y[,,Nlayer] * AuxE2
Epiz.E1[,,Nlayer] <- 0i
## loop downwards to compute all field amplitudes
for (ii in seq(Nlayer-1, 1, by=-1)){
Hiy.H1y[,,ii] <- Mi11[,,ii]*Hiy.H1y[,,ii+1] + Mi12[,,ii]*Hpiy.H1y[,,ii+1]
Hpiy.H1y[,,ii] <- Mi21[,,ii]*Hiy.H1y[,,ii+1] + Mi22[,,ii]*Hpiy.H1y[,,ii+1]
AuxE1 <- matrix(sqrt(epsilon[,1] + 0i) / k0 / epsilon[,ii] , nrow=Nlambda, ncol=Nq)
Eix.E1[,,ii] <- Hiy.H1y[,,ii] * kiz[,,ii] * AuxE1
Epix.E1[,,ii] <- - Hpiy.H1y[,,ii] * kiz[,,ii] * AuxE1
AuxE2 <- outer(epsilon[,1] / epsilon[,ii], Re(q))
Eiz.E1[,,ii] <- - Hiy.H1y[,,ii] * AuxE2
Epiz.E1[,,ii] <- - Hpiy.H1y[,,ii] * AuxE2
}
## loop to compute the local field EFs
for (ii in seq(1, Nlayer-1, by=1)){
## left of interface ii
## the relative coordinate is sampling[[ii]] - thickness[ii]
d1 <- thickness[ii] - sampling[[ii]]
Ml.perp[[ii]] <- sapply(d1, function(.d)
Mod(Eiz.E1[,,ii] * exp(1i*.d*kiz[,,ii]) +
Epiz.E1[,,ii] * exp(-1i*.d*kiz[,,ii]))^2,
simplify="array")
Ml.par[[ii]] <- sapply(d1, function(.d)
Mod(Eix.E1[,,ii] * exp(1i*.d*kiz[,,ii]) +
Epix.E1[,,ii] * exp(-1i*.d*kiz[,,ii]))^2,
simplify="array")
## right of interface ii
## the relative coordinate is sampling[[ii+1]]
Mr.perp[[ii]] <- sapply(sampling[[ii+1]], function(.d)
Mod(Eiz.E1[,,ii+1] * exp( 1i*.d*kiz[,,ii+1]) +
Epiz.E1[,,ii+1] * exp(-1i*.d*kiz[,,ii+1]))^2,
simplify="array")
Mr.par[[ii]] <- sapply(sampling[[ii+1]], function(.d)
Mod(Eix.E1[,,ii+1] * exp( 1i*.d*kiz[,,ii+1]) +
Epix.E1[,,ii+1] * exp(-1i*.d*kiz[,,ii+1]))^2,
simplify="array")
## absolute positions on the right of interface ii
distance[[ii+1]] <- sampling[[ii+1]] + interfaces[ii]
} # end loop
fields <- list(Eix.E1=Eix.E1, Epix.E1=Epix.E1,
Eiz.E1=Eiz.E1, Epiz.E1=Epiz.E1)
}
if(polarisation =="s"){
## ############### ##
## TE polarisation ##
## ############### ##
## field amplitudes p is for prime, reference is field in region 1
Eiy.E1y <- Epiy.E1y <- Hix.H1 <- Hpix.H1 <- Hiz.H1 <- Hpiz.H1 <-
array(0+0i, dim=c(Nlambda, Nq, Nlayer))
Eiy.E1y[,,Nlayer] <- transmission
Epiy.E1y[,,Nlayer] <- 0i
AuxH1 <- matrix(1 / (k0 * sqrt(epsilon[,1] + 0i)), nrow=Nlambda, ncol=Nq)
Hix.H1[,,Nlayer] <- -Eiy.E1y[,,Nlayer] * kiz[,,Nlayer] * AuxH1
Hpix.H1[,,Nlayer] <- 0i
AuxH2 <- matrix(Re(q), nrow=Nlambda, ncol=Nq)
Hiz.H1[,,Nlayer] <- Eiy.E1y[,,Nlayer] * AuxH2
Hpiz.H1[,,Nlayer] <- 0i
## loop downwards to compute all field amplitudes
for (ii in seq(Nlayer-1, 1, by=-1)){
Eiy.E1y[,,ii] <- Mi11[,,ii]*Eiy.E1y[,,ii+1] + Mi12[,,ii]*Epiy.E1y[,,ii+1]
Epiy.E1y[,,ii] <- Mi21[,,ii]*Eiy.E1y[,,ii+1] + Mi22[,,ii]*Epiy.E1y[,,ii+1]
Hix.H1[,,ii] <- - Eiy.E1y[,,ii] * kiz[,,ii] * AuxH1
Hpix.H1[,,ii] <- Epiy.E1y[,,ii] * kiz[,,ii] * AuxH1
Hiz.H1[,,ii] <- Eiy.E1y[,,ii] * AuxH2
Hpiz.H1[,,ii] <- Epiy.E1y[,,ii] * AuxH2
}
## loop to compute the local field EFs
for (ii in seq(1, Nlayer-1, by=1)){
## left of interface ii
## the relative coordinate is sampling[[ii]] - thickness[ii]
d1 <- thickness[ii] - sampling[[ii]]
Ml.par[[ii]] <- sapply(d1, function(.d)
Mod(Eiy.E1y[,,ii] * exp(1i*.d*kiz[,,ii]) +
Epiy.E1y[,,ii] * exp(-1i*.d*kiz[,,ii]))^2,
simplify="array")
## right of interface ii
## the relative coordinate is sampling[[ii+1]]
Mr.par[[ii]] <- sapply(sampling[[ii+1]], function(.d)
Mod(Eiy.E1y[,,ii+1] * exp( 1i*.d*kiz[,,ii+1]) +
Epiy.E1y[,,ii+1] * exp(-1i*.d*kiz[,,ii+1]))^2,
simplify="array")
## Ml.perp[[ii]] <- 0 * Ml.par[[ii]] # default value
## Mr.perp[[ii]] <- 0 * Mr.par[[ii]] # default value
## absolute positions on the right of interface ii
distance[[ii+1]] <- sampling[[ii+1]] + interfaces[ii]
} # end loop
fields <- list(Eiy.E1y=Eiy.E1y, Epiy.E1y=Epiy.E1y)
} # end swich polarisation
## T is nt*cos(Ot)*|Et|^2 / ni*cos(Oi)*|Ei|^2
## for s-pol, |Et|^2 / |Ei|^2 = |ts|^2, hence T = nt/ni * cos(Ot)/cos(Oi) * |ts|^2
## for p-pol, |Et|^2 / |Ei|^2 = (ni/nt)^2 * |tp|^2, hence T = ni/nt * cos(Ot)/cos(Oi) * |tp|^2
# ratio of refractive indices
index.ratio <- matrix(Re(sqrt(epsilon[,1])/sqrt(epsilon[,Nlayer])), nrow=Nlambda, ncol=Nq)
# ratio of cosines
qq <- matrix(q, nrow=Nlambda, ncol=Nq, byrow=TRUE)
m <- Re(sqrt(1 - (index.ratio * qq)^2 + 0i)/sqrt(1 - qq^2 + 0i))
if(polarisation == "p"){
rho <- index.ratio
reflection <- -reflection # sign convention was for H
} else {
rho <- 1 / index.ratio
}
dim(transmission) <- dim(m) # case 1-dims were dropped
R <- Mod(reflection)^2
T <- rho * m * Mod(transmission)^2
T <- drop(T)
transmission <- drop(transmission)
## results
list(wavelength=wavelength, k0 = k0,
angle=angle, q=q,
reflection=reflection, transmission=transmission,
R=R, T=T, A = 1 - R - T,
dist=distance, fields = fields,
Ml.perp=lapply(Ml.perp, drop), Ml.par=lapply(Ml.par, drop),
Mr.perp=lapply(Mr.perp, drop), Mr.par=lapply(Mr.par, drop))
}
##' Multilayer Fresnel coefficients
##'
##' solves the EM problem of a multilayered interface
##' @title multilayercpp
##' @export
##' @param wavelength [vector] wavelength in nm
##' @param k0 [vector] wavevector in nm^-1
##' @param angle [vector] incident angles in radians
##' @param q [vector] normalised incident in-plane wavevector
##' @param epsilon list of N+2 dielectric functions, each of length 1 or length(wavelength)
##' @param thickness vector of N+2 layer thicknesses, first and last are dummy
##' @param ... unused
##' @return fresnel coefficients and field profiles
##' @author baptiste Auguie
##' @examples
##' library(planar)
##' demo(package="planar")
multilayercpp <- function(wavelength = 2*pi/k0, k0 = 2*pi/wavelength,
angle = asin(q), q = sin(angle),
epsilon = list(incident=1.5^2, 1.33),
thickness = c(0, 0), ...){
kx <- outer(k0*sqrt(epsilon[[1]]), q) # kx = q*k0
epsilon = do.call(cbind, epsilon)
## case pure scalars
if(nrow(epsilon) == 1L)
epsilon <- matrix(epsilon, nrow=length(k0),
ncol=length(thickness), byrow=TRUE)
Nlayer <- length(thickness)
Nlambda <- length(k0)
Nq <- length(q)
## checks
stopifnot(thickness[1]==0L,
thickness[Nlayer]==0L)
stopifnot(Nlayer == ncol(epsilon),
nrow(epsilon) == length(k0),
nrow(kx) == length(k0),
ncol(kx) == length(q))
## call the C++ function
res <- cpp_multilayer(as.vector(k0), as.matrix(kx), as.matrix(epsilon),
as.vector(thickness),
0.0, # z=0 irrelevant don't calculate fields
0.0, # psi=0 irrelevant don't calculate fields
FALSE) # don't need intensities
ts <- drop(res$ts)
rs <- drop(res$rs)
tp <- drop(res$tp)
rp <- drop(res$rp)
## T is nt*cos(Ot)*|Et|^2 / ni*cos(Oi)*|Ei|^2
## for s-pol, |Et|^2 / |Ei|^2 = |ts|^2, hence T = nt/ni * cos(Ot)/cos(Oi) * |ts|^2
## for p-pol, |Et|^2 / |Ei|^2 = (ni/nt)^2 * |tp|^2, hence T = ni/nt * cos(Ot)/cos(Oi) * |tp|^2
# ratio of refractive indices
index.ratio <- matrix(Re(sqrt(epsilon[,1])/sqrt(epsilon[,Nlayer])),
nrow=Nlambda, ncol=Nq)
# ratio of cosines
qq <- matrix(q, nrow=Nlambda, ncol=Nq, byrow=TRUE)
m <- Re(sqrt(1 - (index.ratio * qq)^2 + 0i)/sqrt(1 - qq^2 + 0i))
dim(ts) <- dim(tp) <- dim(m) # case 1-dims were dropped
# p-pol
rhop <- index.ratio
rp <- -rp # sign convention was for H
# s-pol
rhos <- 1 / index.ratio
Rp <- Mod(rp)^2
Rs <- Mod(rs)^2
Ts <- drop(rhos * m * Mod(ts)^2)
Tp <- drop(rhop * m * Mod(tp)^2)
ts <- drop(ts)
tp <- drop(tp)
As <- 1 - Rs - Ts
Ap <- 1 - Rp - Tp
list(wavelength=wavelength, k0 = k0,
angle=angle, q=q,
rs=rs, ts=ts,
rp=rp, tp=tp,
Rs=Rs, Ts=Ts, As=As,
Rp=Rp, Tp=Tp, Ap=Ap)
}
##' Multilayer Fresnel coefficients
##'
##' solves the EM problem of a multilayered interface
##' @title multilayerfull
##' @export
##' @param wavelength [vector] wavelength in nm
##' @param k0 [vector] wavevector in nm^-1
##' @param angle [vector] incident angles in radians
##' @param q [vector] normalised incident in-plane wavevector
##' @param epsilon list of N+2 dielectric functions, each of length 1 or length(wavelength)
##' @param thickness vector of N+2 layer thicknesses, first and last are dummy
##' @param psi [numeric] polarisation angle
##' @param z [vector] positions to calculate the electric field intensity
##' @param ... unused
##' @return fresnel coefficients and field profiles
##' @author baptiste Auguie
multilayerfull <- function(wavelength = 2*pi/k0, k0 = 2*pi/wavelength,
angle = asin(q), q = sin(angle),
epsilon = list(incident=1.5^2, 1.33),
thickness = c(0, 0),
psi = 0, z=0, ...){
kx <- outer(k0*sqrt(epsilon[[1]]), q) # kx = q*k0
epsilon = do.call(cbind, epsilon)
Nlayer <- length(thickness)
Nlambda <- length(k0)
Nq <- length(q)
## checks
stopifnot(thickness[1]==0L,
thickness[Nlayer]==0L)
stopifnot(Nlayer == ncol(epsilon),
nrow(epsilon) == length(k0),
nrow(kx) == length(k0),
ncol(kx) == length(q))
## call the C++ function
res <- cpp_multilayer(as.vector(k0),
as.matrix(kx),
as.matrix(epsilon),
as.vector(thickness),
as.vector(z),
as.double(psi), TRUE)
ts <- drop(res$ts)
rs <- drop(res$rs)
tp <- drop(res$tp)
rp <- drop(res$rp)
I <- drop(res$I)
## T is nt*cos(Ot)*|Et|^2 / ni*cos(Oi)*|Ei|^2
## for s-pol, |Et|^2 / |Ei|^2 = |ts|^2,
## hence T = nt/ni * cos(Ot)/cos(Oi) * |ts|^2
## for p-pol, |Et|^2 / |Ei|^2 = (ni/nt)^2 * |tp|^2,
## hence T = ni/nt * cos(Ot)/cos(Oi) * |tp|^2
# ratio of refractive indices
ratio <- matrix(Re(sqrt(epsilon[,1])/sqrt(epsilon[,Nlayer])),
nrow=Nlambda, ncol=Nq)
# ratio of cosines
qq <- matrix(q, nrow=Nlambda, ncol=Nq, byrow=TRUE)
m <- Re(sqrt(1 - (ratio * qq)^2 + 0i)/sqrt(1 - qq^2 + 0i))
rhop <- sqrt(ratio)
rhos <- 1 / rhop
dim(ts) <- dim(tp) <- dim(m) # case 1-dims were dropped
R <- Mod(cos(psi)*rp + sin(psi)*rs)^2
T <- m * Mod(cos(psi)*rhop*tp + sin(psi)*rhos*ts)^2
T <- drop(T)
tp <- drop(tp)
ts <- drop(ts)
A <- 1 - R - T
list(wavelength=wavelength, k0 = k0,
angle=angle, q=q,
rs=rs, ts=ts,
rp=rp, tp=tp,
R=R, T=T, A=A,
I = I)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.