update web

nano-optics · Jan 31, 2022 · a58961f · a58961f
1 parent 7ec0cff
commit a58961f
Show file tree

Hide file tree

Showing 78 changed files with 4,169 additions and 1,480 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/CITATION.cff b/CITATION.cff
@@ -0,0 +1,13 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+title: "mie"
+abstract: Mie scattering
+authors:
+- family-names: "Auguié"
+  given-names: "Baptiste"
+  orcid: "https://orcid.org/0000-0002-2749-5715"
+version: 0.1.0
+date-released: 2021-12-30
+url: "http://nano-optics.ac.nz/mie"
+license: MPL-2.0
+repository-code: "https://github.com/nano-optics/mie/"
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -4,12 +4,12 @@ Title: Mie scattering
 Version: 1.0
 Date: 2014-06-22
 Authors@R: c(person("Baptiste", "Auguie",
-    email = "baptiste.auguie@gmail.com",
+    email = "baptiste.auguie@vuw.ac.nz",
     role = c("aut", "cre"), comment = c(ORCID="0000-0002-2749-5715")))
-License: GPL-3
+License: Mozilla Public License Version 2.0
 Description: Numerical implementation of Mie scattering theory for light
     scattering by spherical particles.
-URL: https://github.com/plasmonics/mie
+URL: https://github.com/nano-optics/mie
 DOI: http://dx.doi.org/10.5281/zenodo.11421
 VignetteBuilder: knitr
 Depends:
@@ -24,4 +24,4 @@ Suggests:
     knitr,
     rmarkdown
 LazyLoad: yes
-RoxygenNote: 5.0.1.9000
+RoxygenNote: 7.1.2
diff --git a/LICENSE b/LICENSE
diff --git a/NAMESPACE b/NAMESPACE
@@ -3,8 +3,11 @@
 export(average_Mloc)
 export(efficiencies)
 export(mie)
+export(mie_approximation)
 export(mie_bh)
+export(mie_ml)
 export(psi)
+export(susceptibilities)
 export(susceptibility)
 export(xi)
 import(Bessel)

diff --git a/R/approximation.R b/R/approximation.R
@@ -15,8 +15,6 @@
 ##' @author Baptiste Auguie
 ##' @family user
 ##' @export
-##' @examples
-##' mie_approximation(50)
 mie_approximation <- function(radius, wavelength, epsilon, medium=1.33, order=Inf, efficiency = FALSE, ...){
 
   s <- sqrt(epsilon) / medium

diff --git a/R/mie-splac.r b/R/mie-splac.r
@@ -1,4 +1,151 @@
 
+##' 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 nmax 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)
+##'
+##'gold <- epsAu(seq(200, 1500))
+##'library(ggplot2)
+##'
+##'params <- expand.grid(order = c(1, 2, Inf), mode = c("EM", "Magnetic", "Electric"), stringsAsFactors=FALSE)
+##'
+##'all <- plyr::mdply(params, mie, wavelength=gold$wavelength, 
+##'              epsilon=gold$epsilon, radius=80, medium=1.5,
+##'             .progress="text")
+##'
+##'m <- reshape2::melt(all, meas = c("extinction", "scattering", "absorption"))
+##'
+##'ggplot(m) +
+##'  facet_grid(mode~variable, 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")
+##' 
+mie <- function(wavelength, epsilon, radius, medium = 1.0,
+                nmax=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 nmax evaluated now
+  coeffs <- susceptibility(s, x, nmax)
+  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)
+}
+
+
+##' Average Mloc
+##'
+##' Enhancement factor averaged over the surface 
+##' @title average_Mloc
+##' @param wavelength real vector
+##' @param epsilon complex vector
+##' @param radius scalar
+##' @param medium scalar, refractive index of surrounding medium
+##' @param nmax truncation order
+##' @return data.frame with wavelength and Mloc
+##' @author Baptiste Auguie
+##' @export
+average_Mloc <- function(wavelength, epsilon, radius, medium = 1.0, nmax=10){
+
+  s <- sqrt(epsilon) / medium
+  x <- 2 * pi / wavelength * medium * radius
+
+  rbx <- ricatti_bessel(x, nmax)
+  GD <- susceptibility(s, x, nmax)
+
+  tmp1 <- rbx[["psi"]] + GD[["G"]] * rbx[["xi"]]
+  tmp2 <- rbx[["dpsi"]] + GD[["D"]] * rbx[["dxi"]]
+  summat1 <- abs(tmp1)^2 + abs(tmp2)^2
+
+  nlen <- seq_len(nmax)
+  cc1 <- t(2*nlen + 1)
+  cc2 <- cc1*nlen*(nlen+1)
+
+  tmp1 <- rbx[["psi"]] + GD[["D"]] * rbx[["xi"]] 
+  summat2 <- abs(tmp1)^2
+
+  # From Eq. H.79
+  data.frame(wavelength=wavelength, Mloc = 1/(2*x^2) * tcrossprod(summat1, cc1) + 1/(2*x^4) * tcrossprod(summat2, cc2))
+}
+
+
+##' 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)
+
+  nmax <- NCOL(GD$G)
+  nvec <- seq_len(nmax)
+  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)
+}
+
 
 ##' Riccati-Bessel function psi and its derivative
 ##'
@@ -96,150 +243,3 @@ susceptibility <- function(s, x, nmax){
        A = 1i * smat   / A_denominator,
        B = 1i * smat   / B_denominator)
 }
-
-##' Average Mloc
-##'
-##' Enhancement factor averaged over the surface 
-##' @title average_Mloc
-##' @param wavelength real vector
-##' @param epsilon complex vector
-##' @param radius scalar
-##' @param medium scalar, refractive index of surrounding medium
-##' @param nmax truncation order
-##' @return data.frame with wavelength and Mloc
-##' @author Baptiste Auguie
-##' @export
-average_Mloc <- function(wavelength, epsilon, radius, medium = 1.0, nmax=10){
-
-  s <- sqrt(epsilon) / medium
-  x <- 2 * pi / wavelength * medium * radius
-
-  rbx <- ricatti_bessel(x, nmax)
-  GD <- susceptibility(s, x, nmax)
-
-  tmp1 <- rbx[["psi"]] + GD[["G"]] * rbx[["xi"]]
-  tmp2 <- rbx[["dpsi"]] + GD[["D"]] * rbx[["dxi"]]
-  summat1 <- abs(tmp1)^2 + abs(tmp2)^2
-
-  nlen <- seq_len(nmax)
-  cc1 <- t(2*nlen + 1)
-  cc2 <- cc1*nlen*(nlen+1)
-
-  tmp1 <- rbx[["psi"]] + GD[["D"]] * rbx[["xi"]] 
-  summat2 <- abs(tmp1)^2
-
-  # From Eq. H.79
- data.frame(wavelength=wavelength, Mloc = 1/(2*x^2) * tcrossprod(summat1, cc1) + 1/(2*x^4) * tcrossprod(summat2, cc2))
-}
-
-
-##' 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)
-
-  nmax <- NCOL(GD$G)
-  nvec <- seq_len(nmax)
-  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 nmax 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)
-##'
-##'gold <- epsAu(seq(200, 1500))
-##'library(ggplot2)
-##'
-##'params <- expand.grid(order = c(1, 2, Inf), mode = c("EM", "Magnetic", "Electric"), stringsAsFactors=FALSE)
-##'
-##'all <- plyr::mdply(params, mie, wavelength=gold$wavelength, 
-##'              epsilon=gold$epsilon, radius=80, medium=1.5,
-##'             .progress="text")
-##'
-##'m <- reshape2::melt(all, meas = c("extinction", "scattering", "absorption"))
-##'
-##'ggplot(m) +
-##'  facet_grid(mode~variable, 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")
-##' 
-mie <- function(wavelength, epsilon, radius, medium = 1.0,
-                nmax=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 nmax evaluated now
-  coeffs <- susceptibility(s, x, nmax)
-  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)
-}
-