inst/tmp/mie.R

In nano-optics/mie: Mie scattering

library(Bessel)

ricatti_bessel <- function(rho, n_max){
  
  rho <- as.complex(rho)
  sq <- sqrt((pi/2)*rho)
  bj <- sq * BesselJ(z = rho, nu = 0.5, nSeq = n_max+1)
  bh <- sq * BesselH(z = rho, nu = 0.5, nSeq = n_max+1, m = 1)
  
  psi <- bj[, -1L]
  xi <- bh[, -1L]
  norho <- outer(1/rho, seq_len(n_max))
  
  dpsi <- bj[, -(n_max+1)] - norho * psi
  dxi <- bh[, -(n_max+1)] - norho * xi
  
  list(psi=psi, xi=xi, dpsi=dpsi, dxi=dxi)
}

rho <- 2*pi/seq(500, 800, by=100) * 50 * sqrt(-10+1i)
n_max <- 10
ricatti_bessel(rho, n_max)

susceptibility <- function(s, x, n_max){
  
  smat <- matrix(s, ncol=n_max, nrow=length(x), byrow=FALSE)
  z <- s*x
  rbz <- ricatti_bessel(z, n_max)
  rbx <- ricatti_bessel(x, n_max)
  PP1 <- rbz[["psi"]] * rbx[["dpsi"]];
  PP2 <- rbx[["psi"]] * rbz[["dpsi"]];
  PP3 <- rbz[["psi"]] * rbx[["dxi"]];
  PP4 <- rbx[["xi"]] * rbz[["dpsi"]];
  
  G_numerator   <-  - PP1 + smat * PP2
  D_numerator   <-    PP2 - smat * PP1
  B_denominator <-  - PP4 + smat * PP3
  A_denominator <-    PP3 - smat * PP4
  
  list(G = G_numerator / A_denominator,
       D = D_numerator / B_denominator,
       A = 1i * smat   / A_denominator,
       B = 1i * smat   / B_denominator)
}


# 
# x <- 2*pi/seq(500, 800, by=100)
# s <- 1.5
# susceptibility1(s, x, 10)
# susceptibility(10, s, x)

##' Efficiencies
##'
##' Calculates the far-field efficiencies for plane-wave illumination
##' @title efficiencies
##' @param x real vector, size parameter
##' @param GD list with Gamma, Delta, A, B
##' @param mode type of mode
##' @param order order of multipoles
##' @return matrix of Qext, Qsca, Qabs
##' @author Baptiste Auguie
##' @export
efficiencies <- function(x, GD, mode=c("EM", "Magnetic", "Electric"), order = NULL){
  
  mode <- match.arg(mode)
  
  n_max <- NCOL(GD$G)
  nvec <- seq.int(n_max)
  nvec2 <- 2 * nvec + 1
  if(all(is.numeric(order)) && all(is.finite(order))) {
    nvec2[-order] <- 0
  }
  
  G2 <- Mod(GD$G)^2
  D2 <- Mod(GD$D)^2
  
  GR <- Re(GD$G)
  DR <- Re(GD$D)
  
  if(mode == "EM"){
    scatmat <- G2 + D2
    ext_coeff <- GR + DR
  } else {
    if(mode == "Electric"){
      scatmat <- D2
      ext_coeff <- DR
    } else {
      scatmat <- G2 
      ext_coeff <- GR
    }
  }
  
  Qsca <- 2 / x^2 * scatmat %*% nvec2
  Qext <- - 2 / x^2 * ext_coeff %*% nvec2
  Qabs <- Qext - Qsca
  
  cbind(Qext = Qext, Qsca = Qsca, Qabs = Qabs)
}

##' Far-field cross-sections
##'
##' Homogeneous sphere illuminated by a plane wave
##' @title mie
##' @param wavelength real vector
##' @param epsilon complex vector
##' @param radius scalar
##' @param medium scalar, refractive index of surrounding medium
##' @param n_max truncation order
##' @param efficiency logical, scale by geometrical cross-sections
##' @param mode type of mode
##' @param order order of multipoles
##' @return data.frame
##' @author Baptiste Auguie
##' @family user
##' @export
##' @examples 
##' gold <- epsAu(seq(400, 800))
##' cross_sections <- with(gold, mie(wavelength, epsilon, radius=50, medium=1.33, efficiency=FALSE))
##' matplot(cross_sections$wavelength, cross_sections[, -1], type="l", lty=1,
##'         xlab=expression(lambda/mu*m), ylab=expression(sigma/mu*m^2))
##' legend("topright", names(cross_sections)[-1], col=1:3, lty=1)
mie <- function(wavelength, epsilon, radius, medium = 1.0,
                n_max=ceiling(2 + max(x) + 4 * max(x)^(1/3)),
                efficiency = FALSE, mode=c("EM", "Magnetic", "Electric"),
                order = Inf){
  
  mode <- match.arg(mode)
  
  s <- sqrt(epsilon) / medium
  x <- 2 * pi / wavelength * medium * radius
  
  ## lazy evaluation rules.. default n_max evaluated now
  coeffs <- susceptibility(s, x, n_max)
  Q <- efficiencies(x, coeffs, mode=mode, order=order)
  if(!efficiency) Q <- Q * (pi*radius^2)
  results <- data.frame(wavelength, Q)
  names(results) <- c("wavelength", "extinction", "scattering", "absorption")
  invisible(results)
}

library(dielectric)
gold <- epsAg(seq(300, 800))
a <- 30
cross_sections <- with(gold, mie(wavelength, epsilon, radius=a, medium=1.33, efficiency=TRUE, order=1))
matplot(cross_sections$wavelength, cross_sections[, -1], type="l", lty=1,
        xlab=expression(lambda/mu*m), ylab=expression(sigma/mu*m^2))
# M <- with(gold, average_Mloc(wavelength, epsilon, radius=a, medium=1.33, n_max=100))
# lines(M$wavelength, M$Mloc/max(M$Mloc)*max(cross_sections[,-1]), lty=2)
legend("topright", c(names(cross_sections)[-1], "<Mloc>"), col=1:3, lty=c(1,1,1,2))

 
gold <- epsAu(seq(200, 1500))
library(ggplot2)

params <- expand.grid(order = c(1, 2, 3, Inf), mode = c("EM", "Magnetic", "Electric"), stringsAsFactors=FALSE)

all <- purrr::pmap_df(params, mie, wavelength=gold$wavelength, 
                   epsilon=gold$epsilon, radius=80, medium=1.5,
                   .id = 'id')

params$id <- as.character(1:nrow(params))

m <- tidyr::pivot_longer(left_join(params, all, by='id'),
                         cols = c("extinction", "scattering", "absorption"))

ggplot(m) +
  facet_grid(mode~name, scales="free") +
  geom_path(aes(wavelength, value, colour = factor(order))) +
  scale_linetype_manual(values = c(2, 3, 1)) +
  labs(x = expression(wavelength / nm),
       y = expression(sigma / nm^2),
       colour = "Mode",
       linetype = "Order")