## Surface plasmon-polaritons in the Kretschmann configuration
_baptiste Auguié -- `r format(Sys.time(), "%d %B, %Y")`_

The Kretschmann configuration consists of a thin metal film, typically 50nm of gold or silver, deposited on top of a high-index prism (n=1.5 for glass). Light incident from the prism side undergoes total internal reflection (TIR) above ~45 degrees (internal angle). The evanescent field associated with TIR penetrates the metal and may couple to surface plasmon-polaritons supported at the air/metal interface.

Here we model the optical properties of such a system, starting with the angular variation of the reflectivity.

### Reflectivity against internal incident angle for the Kretschmann configuration, at fixed wavelength

```{r , echo=FALSE,results='hide'}
library(knitr)
library(ggplot2)
opts_chunk$set(fig.path="kretschmannreflectivity/", fig.width=10,
               warning=FALSE,error=FALSE,message=FALSE,tidy=TRUE)
library(ggplot2)
theme_set(theme_minimal() + theme(panel.border=element_rect(fill=NA)))
library(RColorBrewer)
col <- brewer.pal(3,"PRGn")
```

#### Setting up

```{r , results='hide'}
library(planar)
library(ggplot2)
require(reshape2)
library(gridExtra)
require(plyr)
```

#### Modelling the reflectivity

```{r}
wvl <- 632.8
gold <- epsAu(wvl)
results <- recursive_fresnelcpp(epsilon=list(1.5^2, gold$epsilon, 1.0),
                                wavelength=gold$wavelength, thickness=c(0, 50, 0),
                                angle=seq(0, pi/2, length=2e3), polarisation='p')
str(results)
```

#### Plotting the results

```{r reflectivity,fig.width=10}
m <- data.frame(results[c("angle", "R")])

tir <- asin(1/1.5) * 180/pi
 
ggplot(m) +
  geom_vline(aes(xintercept=x),
             data=data.frame(x=tir),
             linetype=2,color="grey50") +
  geom_line(aes(angle*180/pi, R)) +
  scale_y_continuous("Reflectivity", expand=c(0,0), limits=c(0,1))+
  scale_x_continuous("Internal angle /degrees", expand=c(0,0), 
                     breaks=seq(0,90, by=15)) 
  
```

#### Variation of the parameters, and effect on the resonance

We now look at the effect of changing the thickness of the metal layer, from non-existent (single air/glass interface), to an opaque metal film. First, we wrap the calculation in a function, and loop over this function with a vector of film thicknesses.

```{r loop}

simulation <- function(thickness = 50){
results <- recursive_fresnelcpp(epsilon=list(1.5^2, gold$epsilon, 1.0^2),
                                wavelength=gold$wavelength, 
                                thickness=c(0, thickness, 0),
                                angle=pi/180*seq(15, 60, length=500), 
                                polarisation='p')
data.frame(results[c("angle", "R")])

}

## loop over parameters
parameters <- function(res=10) 
  data.frame(thickness = seq(0, 100, length=res))

d1 <- mdply(parameters(10), simulation)
d2 <- mdply(parameters(300), simulation)


p1 <- 
ggplot(d1) +
  geom_line(aes(angle*180/pi, R, colour=thickness, group=thickness)) +
  scale_y_continuous("Reflectivity", expand=c(0,0), limits=c(0,1))+
  scale_x_continuous("Internal angle /degrees", expand=c(0,0), 
                     breaks=seq(0,90, by=15)) +
  guides(colour=guide_legend()) 

## colour map
p2 <- 
ggplot(d2) +
  geom_raster(aes(angle*180/pi, thickness, fill=R)) +
  scale_y_continuous("thickness", expand=c(0,0))+
  scale_x_continuous("Internal angle /degrees", expand=c(0,0), 
                     breaks=seq(0,90, by=15)) 

grid.arrange(p1, p2, nrow=2)
```


```{r variation}
minimum <- ddply(subset(d2, angle > tir*pi/180 & thickness > 5), .(thickness), summarize, 
                 angle = angle[which.min(R)] * 180/pi,
                 min = min(R))
ggplot(melt(minimum, id="thickness")) + 
  facet_grid(variable~., scales="free") +
  geom_line(aes(thickness, value)) +
  labs(y="", x="thickness /nm")
  
```
