---
title: "RLumModel - Using own parameter sets"
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 - Using own parameter sets}         
  %\VignetteEncoding{UTF-8}{inputenc}
  %\VignetteEngine{knitr::rmarkdown} 
---
```{r, echo=FALSE, message = FALSE}
library(RLumModel)
```

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.pos = 'H', fig.align = 'center')
```

# Introduction

With 'RLumModel' &ge; 0.2.0 is it possible to simulate quartz luminescence behaviour of own parameters or parameter sets, which are not included in the package but also known in literature. Wide spreaded OTOR (One-Trap-One-Recombination-center) models can be included, too. 
This vignette gives three comprehensive examples how to implement parameter sets and prooves the results recalculating the original simulations.

'RLumModel' offers maximum flexibility and fast calculation of ordinary first-order differential equations (ODEs) describing luminescence behaviour, because of:

* flexible handling of different numbers of electron traps and hole centres in generating the model equations. The user has not to care about coding ODEs in **R**. 
* Solving these equations is performed with C++ [@Eddelbuettel_2013] and offers a fast calculation.

# Examples

This chapter shows the handling of own parameter sets with 'RLumModel'. For this purpose three model parameters known from literature were taken: @Pagonis_2009a, @Lawless_2009 and @Chen_2013d. 

## Pagonis 2009 

@Pagonis_2009a presented parameters for their luminescence modelling of radioluminescence. 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 manual for involving these parameters in 'RLumModel' and re-calculating the simulationa made by @Pagonis_2009a. Note that in the original publication Figure 3 and Figure 6 are inconsistent with each other. For a doserate of 0.1 Gy/s an initial RL intensity of ca. 1.5e14 is obtained (see Figure 6 in original publication and simulations below). 

### Set parameters

First of all the model parameters had to be set. In 'RLumModel' this can be done via \texttt{list()}.
The list has to contain the following items:

* N: Concentration of electron- and hole traps [cm<sup>-3</sup>]
* E: Electron/Hole trap depth [eV]
* s: Frequency factor [s<sup>-1</sup>]
* A: Conduction band to electron trap and valence band to hole trap transition probability [s<sup>-1</sup> cm <sup>3</sup>].
* B: Conduction band to hole centre transition probability [s<sup>-1</sup> cm <sup>3</sup>].
* Th: Photo-eviction constant or photoionisation cross section, respectively. If not set: by default 0.
* E_th: Thermal assistence energy [eV]. If not set: by default 0.
* k_B: Boltzman constant 8.617e-05 [eV/K]. If not set: by default  8.617e-05.
* W: activation energy 0.64 [eV] (for UV). If not set: by default  0.64.
* K: 2.8e7 (dimensionless constant). If not set: by default  2.8e7.
* model: "customized"
* R (optional): Ionisation rate (pair production rate) equivalent to 1 Gy/s [s<sup>-1</sup> cm <sup>-3</sup>]

**Note**: 

* Not every publication uses the same definition of parameter `A` and `B`.
* When no thermal quenching is expected, set `K = 0`. A numerical value of `W` is necessary, if \texttt{K $\neq$ 0}.
* The parameter `model = "customized"` is necessary to not load a stored model within the RLumModel pacakge.
* The luminescence center, which is responsible for the luminescence production (TL-, OSL-, RF-signal), has to be the last entry in the model parameters, see examples below.
* For further details of the parameter see @Bailey_2001 and @Wintle_1975.
* The first two entries in `N`, `A` and `B` belong to the electron traps and the last two entries to the hole centres. This order is necessary.
* The first entry in `N`, `E`, `s`, `A`, `B`, `Th` and `E_th` belong to the first energy level (e.g. 110°C), the second, third, ... entries to second, third, ... energy levels.
* The entries `Th` and `E_th` are not necessary, because when they miss in the definition of the own parameters they will automatically be set to 0. 

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

It is important to notice, that in @Pagonis_2009a `B` is the valence band to hole centre probability, but in @Bailey_2001 this is A<sub>j</sub>. The default setting of RLumModel is the definition by @Bailey_2001 and so the values of `B` (in @Pagonis_2009a) are `A` in the notation above.

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

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

  
### Running the simulation with RLumModel

When all parameters are set, the simulation can be started. The main function in RLumModel is `model_LuminescenceSignals()` and the usage with own parameter sets is described below. For a general overview for creating a sequence, running RLumModel with stored models etc. the user is referred to @Friedrich_2016 and to the vignette [RLumModel - Getting started with RLumModel](RLumModel_-_Getting_started_with_RLumModel.html).

For simulating the results of @Pagonis_2009a the follwing sequence is needed.

```{r set sequence Pagonis 2009}  
sequence <- list(RF = c(20, 0.1, 0.1))
```
  
This sequence describes a radiofluorescence simulation at 20 °C with a dose of 0.1 Gy and a dose rate of 0.1 Gy/s, so the stimulation time is 1s. 

The parameters `own_parameters` and `own_state_parameters` in `model_LuminescenceSignals()` are prepared for using own created parameter sets. Parameter `model =  "customized"` is necessary to not load stored parameters. 
  
  
```{r, fig.cap = "RF signal for 0.1 Gy/s"}  
RF_Pagonis2009 <- model_LuminescenceSignals(
  model = "customized", 
  sequence = sequence, 
  own_parameters = own_parameters, 
  own_state_parameters = own_state_parameters,
  verbose = FALSE)
```
  
As in the original publication, initially the RF signal increases and is followed by an approximately linear region until the stimulation ends.
Figure 5 in @Pagonis_2009a shows the concentration of the luminescence center m1 for the stimulation time of 1s.
With RLumModel this can be plotted very fast with the following command (for a detailed description see vignette [RLumModel - Getting started with RLumModel](RLumModel_-_Getting_started_with_RLumModel.html)
  
```{r, fig.cap = "Concentration of m1 during RF"}
concentration_m1 <- Luminescence::get_RLum(
  RF_Pagonis2009,
  recordType = c("conc. level 4"))
  
Luminescence::plot_RLum(
  concentration_m1, 
  ylim = c(9.2e15, 9.6e15))
```


### Re-calculate the original results

Reproducing Figure 3 and Figure 6 in @Pagonis_2009a a loop over different dose rates is necessary. 
The following code lines are able to run the model for five different dose rates from 0.1 to 0.5 Gy/s and plot all contained RF curves and the initial RF signal.
For a more detailed descripton of the loop and the single commands therein the user is referred to @Friedrich_2016 and the vignette [RLumModel - Getting started with RLumModel](RLumModel_-_Getting_started_with_RLumModel.html).

```{r Different Dose-rates Pagonis 2009, fig.cap = "RF signals for different dose rates"}
dose.rate <- seq(from = 0.1, to = 0.5, by = 0.1)

model.output <- lapply(dose.rate, function(x) {
    
    sequence <- list(RF = c(20, x, x))
    
    RF_data <- model_LuminescenceSignals(
      model = "customized", 
      sequence = sequence, 
      own_parameters = own_parameters, 
      own_state_parameters = own_state_parameters,
      verbose = FALSE,
      plot = FALSE
    )
    
    ## "RF$" for exact matching RF and not (RF)
    return(get_RLum(RF_data, recordType = "RF$", drop = FALSE))
    
  })

model.output.merged <- merge_RLum(model.output)

plot_RLum(
 object = model.output.merged,
 xlab = "Stimulation time [s]",
 ylab = "RF signal [a.u.]",
 legend.text = paste(dose.rate, "Gy/s"),
 legend.pos = "outside",
 combine = TRUE)
```



The following code calcultes the initial RF signal for the five different dose rates.


```{r Initial signal}
dose.rate <- seq(from = 0.1, to = 0.5, by = 0.1)

model.output <- vapply(X = dose.rate, FUN = function(x) {
    
    sequence <- list(RF = c(20, x, x))
    
    temp <- model_LuminescenceSignals(
      model = "customized", 
      sequence = sequence, 
      own_parameters = own_parameters, 
      own_state_parameters = own_state_parameters,
      verbose = FALSE,
      plot = FALSE
    )
    
    ## "RF$" for exact matching RF and not (RF)
    RF_curve <- get_RLum(temp, recordType = "RF$")
    
    return(max(get_RLum(RF_curve)[2,2]))
    
  }, FUN.VALUE = 1)
```

```{r, echo=FALSE, fig.cap = "Initial RF signal for different dose rates with parameters of Lawless 2009"}
plot(
  dose.rate,
  model.output,
  type = "b",
  xlab = "Stimulation Time [s]",
  ylab = "Initial RF intensitiy [a.u.]"
  )
```


The results show that 'RLumModel' is able to simulate the same results as published in @Pagonis_2009a with only little effort. All these examples can be modified to own needs, e.g. own sequences or own parameters.

## Lawless 2009

@Lawless_2009 investigateted the sublinear dose dependence of TL and OSL. They published a set of model parameters to simulate the behaviour of the quartz luminescence
system during different dose rates. In contrast to the example above, this simulation has no state parameters and so they were not definded.

### Set parameters and recalculate the results

All used parameters are defined in the named list `own_parameters`. `K=0` was chosen, because no thermal quenching was simulated.
Note: In the "Bailey 2001" notation `B` has the same meaning as `Am` in @Lawless_2009 (for details see example in chapter 2.1.1). 

```{r Lawless 2009 set parameters}
own_parameters <- list(
  N = c(1e14, 1e15),
  E = c(0, 0),
  s = c(0, 0),
  A = c(1e-13, 1e-14),
  B = c(0, 1e-7),
  K = 0,
  model = "customized",
  R = 1e8)

sequence <- list(RF = c(20, 100, 1))

RF_Lawless_2009 <- model_LuminescenceSignals(
  model = "customized", 
  sequence = sequence, 
  own_parameters = own_parameters,
  verbose = FALSE,
  plot = FALSE)

concentration_n <- Luminescence::get_RLum(
  RF_Lawless_2009, 
    recordType = c("conc. level 1"))
```

This code leads to the following results and shows the same as plotted in @Lawless_2009, Fig. 2 (plot commands not shown here). More details to the equations mentioned in the legend are available in the original publication.

```{r, echo=FALSE, fig.cap = "Concentration of Level 1 with numerical and analytical solutions"}
Luminescence::plot_RLum(
  concentration_n,
  ylim = c(0, 15e8), lwd = 3)

t <- seq(0, 100, 2)
numerical_eq16 <- 1e-13*1e14/1e-7 *((1 + 2*1e-7*1e8*t/(1e-13*1e14))^(0.5)-1)

numerical_eq18 <- (2*1e-13*1e14*1e8*t/(1e-7))^(0.5)

lines(t, numerical_eq16, pch = 3, col = "red", type = "b")
lines(t, numerical_eq18, pch = 4, col = "green", type = "b")

legend("bottomright", legend = c("Simulated", "Eq. 16","Eq. 18"), col = c("black", "red", "green"), lwd = 1)
```




## Chen 2013

@Chen_2013d published a numerical model to investigate the quasi-equilibrium assumptions in TL. For the description of the system a OTOR model was used. 


### Set parameters

This model is the first in this vignette which did not start its simulation at 20 °C. For this cases, `model_LuminescenceSignals()` offers a parameter called `own_start_temperature`. This parameter offers maximal flexibility for the user to set the initial temperature of the simulation. The parameter takes effect when `model = "customized"` is used, see example below.

```{r Chen 2013 set parameters}
own_parameters <- list(
  N = c(1e9, 0),
  E = c(0.4, 0),
  s = c(1e11, 0),
  A = c(1e-9,0),
  B = c(0, 1e-10),
  K = 0,
  model = "customized")

own_state_parameters <- c(1e8, 1e8)

own_start_temperature <- -220

sequence <- list(TL = c(-220, 130, 1))
```


### Re-calculate the original results

Here the parameter `own_start_temperature` from the function `model_LuminescenceSignals()` is used to set the beginning of the measurement to -220°C. It is important, that 'RLumModel' always uses temperatures in °C.

```{r solve Chen 2013, fig.cap = "TL with parameter sets of Chen 2013"}
TL_Chen2013 <- model_LuminescenceSignals(
  model = "customized", 
  sequence = sequence, 
  own_parameters = own_parameters, 
  own_state_parameters = own_state_parameters,
  own_start_temperature = own_start_temperature,
  verbose = FALSE)
```

With this result it is possible to plot the concentration of every single energy level, leading to the following plot (see also Fig. 6 in @Chen_2013d)

```{r, echo=FALSE, fig.cap = "Concentrations of different energy levels"}
concentration <- Luminescence::get_RLum(
  TL_Chen2013, 
  recordType = c("conc. level 1", "conc. level 2", "conc. n_c"),
  drop = FALSE)

concentration@records[[1]]@recordType <- "TL"
concentration@records[[2]]@recordType <- "TL"
concentration@records[[3]]@recordType <- "TL"

Luminescence::plot_RLum(
  concentration,
  combine = TRUE,
  ylab = "concentrations",
  main = "",
  legend.text = c("n", "m","nc")
  )
```


# Summary

This vignette showed the potential of the R package 'RLumModel' to use own parameter sets simulating quartz luminescence behaviour.
Quartz as well as Al<sub>2</sub>O<sub>3</sub> luminescence phenomena can be numerically described and graphically plotted.

# References