chi <- function(laxis = 50, saxis = 50) {
    e <- sqrt(1 - saxis^2/laxis^2)
    (1 - e^2)/e^2 * (-1 + 1/(2 * e) * log((1 + e)/(1 - e)))
}

polarizability_ellipsoid <- function(wavelength, epsilon, a = 50, b = 30, medium = 1.33) {
    La <- if (a == b) 
        1/3 else chi(a, b)
    Lb <- Lc <- 1/2 * (1 - La)

    Vcm <- a * b * b
    V <- 4 * pi/3 * Vcm

    Vcm * (epsilon - medium^2) * 1/cbind(medium^2 + La * epsilon, medium^2 + 
        Lb * epsilon, medium^2 + Lc * epsilon)
}

extinction <- function(polarizability, k) {
    4 * pi * k * Im(polarizability)
}

drude <- function(wavelength, epsilon.inf = 5, lambda.p = 300, mu.p = 5000) {

    data.frame(wavelength = wavelength, epsilon = epsilon.inf * (1 - 1/(lambda.p^2 * 
        (1/wavelength^2 + (0+1i)/(mu.p * wavelength)))))
}

wavelength <- seq(200, 800, length = 1000)
medium <- 1
metal <- drude(wavelength)
k <- 2 * pi/(medium * wavelength)
a0 <- 20
ar <- 4
# 4/3*pi*a*b^2 = 4/3*pi*a0^3 a = b*ar a = a0 * ar^2/3 b = a / ar = a0 *
# ar^(-1/3)

alpha_sphere <- polarizability_ellipsoid(wavelength, metal$epsilon, a = a0, 
    b = a0, medium = medium)
alpha_rod <- polarizability_ellipsoid(wavelength, metal$epsilon, a = a0 * ar^(2/3), 
    b = a0 * ar^(-1/3), medium = medium)

extinction_sphere <- rowMeans(extinction(alpha_sphere, k))
extinction_rod <- rowMeans(extinction(alpha_rod, k))

matplot(wavelength, cbind(extinction_sphere, extinction_rod), t = "l", lty = 1, 
    col = 1:2, bty = "n", lwd = 3, xlab = "wavelength / nm", ylab = "", yaxt = "n")

pdf("export.pdf", width = 8, height = 5)
library(grid)
grid.newpage()
# '#E41A1C' '#377EB8' '#4DAF4A'
pushViewport(dataViewport(xData = wavelength, yData = cbind(extinction_sphere, 
    extinction_rod)))
grid.lines(wavelength, extinction_sphere, def = "native", gp = gpar(lwd = 3, 
    col = "#E41A1C"))
grid.lines(wavelength, extinction_rod, def = "native", gp = gpar(lwd = 3, col = "#377EB8"))
upViewport()
dev.off()

## pdf 
##   2

svg("export.svg", width = 8, height = 5)
library(grid)
grid.newpage()
# '#E41A1C' '#377EB8' '#4DAF4A'
pushViewport(dataViewport(xData = wavelength, yData = cbind(extinction_sphere, 
    extinction_rod)))
grid.lines(wavelength, extinction_sphere, def = "native", gp = gpar(lwd = 3, 
    col = "#E41A1C"))
grid.lines(wavelength, extinction_rod, def = "native", gp = gpar(lwd = 3, col = "#377EB8"))
upViewport()
dev.off()

## pdf 
##   2