---
title: "RLumModel - Simulating luminescence in Al<sub>2</sub>O<sub>3</sub>"
author: "Johannes Friedrich, University of Bayreuth"
date: '`r Sys.Date()`'
output:
  rmarkdown::html_vignette:
    fig_caption: yes
    fig_height: 5
    fig_width: 5
    number_sections: yes
    toc: yes
csl: elsevier-harvard.csl
bibliography: RLumModel.bib
vignette: |
  %\VignetteIndexEntry{RLumModel - Simulating luminescence in Al<sub>2</sub>O<sub>3</sub>}         
  %\VignetteEncoding{UTF-8}{inputenc}
  %\VignetteEngine{knitr::rmarkdown} 
---
```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.pos = 'H', fig.align = 'center')
```

# Introduction

This vignette describes the usage of the **R**-package RLumModel for simulating Al<sub>2</sub>O<sub>3</sub> behaviour. In literature a lot of parameters for describing the thermoluminescence (TL) and/or opical stimulated luminescence (OSL) are given. We give same examples for simulating different phenomena of luminescence in Al<sub>2</sub>O<sub>3</sub>. Note that most of the model do not support TL **and** OSL simulations. So please be careful when creating sequences. The models presented in Sec. 2.1 and 2.2 support TL simulations, but no OSL because parameters for `E_th` (Thermal assistence energy) and `Th` (Photo-eviction constant or photoionisation cross section).

# Examples

All examples need the **R** package RLumModel, so let's load it!

```{r, message = FALSE}
library(RLumModel)
```

## Akselrod 1998

@Akselrod_1998 presented parameters for modelling the effect of quenching TL signals. This model was built for Al<sub>2</sub>O<sub>3</sub>, but the rate equations are identical with describing electron movements in quartz. Below is a step-by-step introduction for involving these parameters in 'RLumModel' and re-calculating the simulationa from Figure 9 in @Akselrod_1998.

### Set own starting parameters

For a detailed introduction to own parameter sets in RLumModel, see vignette [RLumModel - Using own parameter sets](RLumModel_-_Using_own_parameter_sets.html).

As a next step it is possible to set own starting-parameters, also called state parameters. In the case of @Akselrod_1998 they submitted initial concentrations of electrons and holes. This can be done via:
  
```{r}  
own_state_parameters <- c(0, 0, 1e11)
```

Here the first entry is the first electron trap, the second entry the second electron trap and the third entry the luminescence center responsible for the TL signal.
The vector `own_state_parameters` needs as much entries as energy levels used in the model.


### Calculate thermal quenching effect

The effect of quenching luminescence signals will be simulated. In RLumModel it is possible to submit the parameter `K` as an indicator for recognising thermal quenching or not. If the parameter is set to `0` thermal quenching will be neglected. Otherwise it wil be calculated.
  
  
```{r, fig.cap = "The effect of thermal quenching for TL signals"}  
W <- c(0, 1.05)

TL_Akselrod_1998 <- lapply(W, function(W){
  
  own_parameters <- list(
    N = c(1e18, 1e17, 1e18),
    E = c(1.25, 0, 0),
    s = c(1e13, 0, 0),
    A = c(1e-16, 1e-16, 2e-17),
    B = c(0, 0, 1e-16),
    K = ifelse(W == 0, 0, 1e12),
    W = W,
    model = "customized",
    R = 1e13)
  
  sequence <- list(
    IRR = c(20, 0.225, 0.01),
    TL = c(20, 300, 1))
  
  model_output <- model_LuminescenceSignals(
    model = "customized", 
    sequence = sequence, 
    own_parameters = own_parameters, 
    own_state_parameters = own_state_parameters,
    verbose = FALSE,
    plot = FALSE)
  
  return(get_RLum(model_output, recordType = "TL$", drop = FALSE))
})

merge_results <- merge_RLum.Analysis(TL_Akselrod_1998)

plot_RLum.Analysis(
  merge_results,
  legend.text = c("Unquenched", "Quenched"),
  combine = T)
```
  

## Pagonis 2007

@Pagonis_2007b published three different parameter sets for Al<sub>2</sub>O<sub>3</sub> chips. Here we analyse 'Chip101' and show how to re-calculate the results presented in the publication.

### Set parameters for Chip101

```{r Pagonis_2007 Chip101}
own_parameters <- list(
  N = c(2e15, 2e15, 2.4e16, 1e17),
  E = c(1.3,0, 0, 0),
  s = c(1e13, 0, 0, 0),
  A = c(2e-8, 2e-9, 4e-9, 1e-8),
  B = c(0, 0, 5e-11, 4e-8),
  K = 1e11,
  W = 1.1,
  model = "customized",
  R = 1.7e15)

own_state_parameters <- c(0, 0, 0, 9.4e15)
```


### Caclualte TL signal dependend on the given dose

We now have to define the sequence for simulating the TL behaviour. For that purpose different doses are given before the TL measurement. To handle this a `lapply` command was used. The result is the combined plot of all simulated TL measurements.

```{r, fig.cap= "TL signal for different given doses"}
dose <- 10^seq(-1, 3.5, 0.5)

Pagonis_2007 <- lapply(dose, function(dose){
  
  sequence <- list(
    IRR = c(20, dose, 1),
    PAUSE = c(20, 60),
    TL = c(20, 250, 1))
  
  model_output <- model_LuminescenceSignals(
    model = "customized", 
    sequence = sequence, 
    own_parameters = own_parameters,
    own_state_parameters = own_state_parameters,
    verbose = FALSE,
    plot = FALSE)

  return(Luminescence::get_RLum(model_output, recordType = "TL", drop = FALSE))
})

merge_results <- Luminescence::merge_RLum.Analysis(Pagonis_2007)

Luminescence::plot_RLum.Analysis(
  merge_results,
  subset =  list(recordType = "TL$"),
  xlim = c(100, 250),
  legend.text = paste0(round(dose, digits = 2), " Gy"),
  combine = T)

```

### Calculate maximum TL signal and concentration of L-centre

Figure 2 in the original publiction by @Pagonis_2007b will be plotted with the following commands. The following code commands are able to calculate the maximum of the TL signal of all curves. 

```{r calculate TL max}
TL_max <- vapply(1:length(Pagonis_2007), function(x){
  
  TL <- get_RLum(get_RLum(Pagonis_2007[[x]], recordType = "TL$"))
  
  return(max(TL[,2]))
  
}, FUN.VALUE = 1)
```


Now the calculation of the concentration:

```{r calculate concentration of m1}

m1_max <- vapply(1:length(Pagonis_2007), function(x){
  
  m1 <- get_RLum(get_RLum(Pagonis_2007[[x]], recordType = "conc. level 4"))
  
  return(m1[1,2])

}, FUN.VALUE = 1)
```

When plotting L-centre concentration vs. dose or TL<sub>max</sub> vs. dose the same results as presented in Fig. 2 in @Pagonis_2007b are reached. 

```{r plot TL max, echo = FALSE, fig.cap = "TL max vs. dose"}
plot(dose, TL_max, log = "xy",
     xlab = "Dose [Gy]", ylab = "TL max [a.u.]",
     ylim = c(1e11, 1e14), pch = 16)
grid()
```

```{r plot concentration m1, echo = FALSE, fig.cap = "Concentration of m1 (L-centre) vs. dose"}
plot(x = dose, y = m1_max, log = "x",
     xlab = "Dose [Gy]", ylab = "Conc. m1 [a.u.]",
     ylim = c(1e14, 1.2e16), pch = 16)
grid()
```

### Concentration at the end of irradiation

```{r}
dose <- seq(1,400, 1)

Fig_5 <- lapply(dose, function(dose){
  
  sequence <- list(
    RF = c(20, dose, 1))
  
  model_output <- model_LuminescenceSignals(
    model = "customized", 
    sequence = sequence, 
    own_parameters = own_parameters,
    own_state_parameters = own_state_parameters,
    verbose = FALSE,
    plot = FALSE)

})
```

```{r}
n1 <- vapply(1:length(Fig_5), function(x){
  
  temp <- get_RLum(get_RLum(Fig_5[[x]], recordType = "conc. level 1"))
  
  return(temp[nrow(temp),2])
  
}, FUN.VALUE = 1)

n2 <- vapply(1:length(Fig_5), function(x){
  
  temp <- get_RLum(get_RLum(Fig_5[[x]], recordType = "conc. level 2"))
  
  return(temp[nrow(temp),2])
  
}, FUN.VALUE = 1)

m1 <- vapply(1:length(Fig_5), function(x){
  
  temp <- get_RLum(get_RLum(Fig_5[[x]], recordType = "conc. level 4"))
  
  return(temp[nrow(temp),2])
  
}, FUN.VALUE = 1)

m2 <- vapply(1:length(Fig_5), function(x){
  
  temp <- get_RLum(get_RLum(Fig_5[[x]], recordType = "conc. level 3"))
  
  return(temp[nrow(temp),2])
  
}, FUN.VALUE = 1)
```


```{r, fig.cap = "Concentration in different traps at the end of the irradiation"}
plot(dose, m2, type = "l", ylim = c(0, 1.5e16), xlim = c(0, 500), xlab= "Dose [Gy]", ylab = "Concentration [a.u.]")
lines(dose, n2, col = "red")
lines(dose, m1, col = "green")
lines(dose, n1, col = "blue")

legend("topright", legend = c("n1", "n2", "m1", "m2"), col = c("blue", "red", "green", "black"), lwd = 1, bty = "n")
grid()
```


# Summary

This vignette showed the potential of the **R** package 'RLumModel' to simulate Al<sub>2</sub>O<sub>3</sub> behaviour. Two different models known from literature were re-compiled in the framework of RLumModel.

# References