R-Lum / RLumModel
1
#' sequence step RF/RL simulation
2
#'
3
#' This function simulates the radioluminescence of quartz in the energy-band-model.
4
#'
5
#' @param temp \code{\link{cnumeric}} (\bold{required}): temperature [deg. C] at which the dose should be applied
6
#'
7
#' @param dose \code{\link{numeric}} (\bold{required}): dose to apply in Gray
8
#'
9
#' @param doesRate \code{\link{numeric}} (\bold{required}): named list with model parameters. Note:
10
#' not every parameter apply to every model, see details for further information
11
#'
12
#' @param n \code{\link{numeric}} or \code{\linkS4class{RLum.Results}} (\bold{required}):
13
#' concentration of electron-/holetraps, valence- and conduction band
14
#' from step before. This is necessary to get the boundary condition for the ODEs.
15
#'
16
#' @param RLumModel_ID \code{\link{numeric}} (optional): A ID-number for the RF-step. This ID
17
#' is pass down to \link{calc_concentrations} so all concentrations had the same ID as the
18
#' sequence step they were calculated from. This ID is identic to the sequence step in "sequence".
19
#'
20
#' @param parms \code{\linkS4class{RLum.Results}} (\bold{required}): The specific model parameters are used to simulate
21
#' numerical quartz luminescence results.
22
#'
23
#' @param \dots further arguments and graphical parameters passed to
24
#' \code{\link{plot.default}}. See details for further information
25
#'
26
#' @return This function returns an RLum.Results object of the RF/RL simulation.
27
#'
28
#' @section Function version: 0.1.3 [2017-11-20]
29
#'
30
#' @author Johannes Friedrich, University of Bayreuth (Germany),
31
#'
32
#' @references
33
#'
34
#' Bailey, R.M., 2001. Towards a general kinetic model for optically and thermally stimulated
35
#' luminescence of quartz. Radiation Measurements 33, 17-45.
36
#'
37
#' Bailey, R.M., 2002. Simulations of variability in the luminescence characteristics of natural
38
#' quartz and its implications for estimates of absorbed dose.
39
#' Radiation Protection Dosimetry 100, 33-38.
40
#'
41
#' Bailey, R.M., 2004. Paper I-simulation of dose absorption in quartz over geological timescales
42
#' and it simplications for the precision and accuracy of optical dating.
43
#' Radiation Measurements 38, 299-310.
44
#'
45
#' Pagonis, V., Chen, R., Wintle, A.G., 2007: Modelling thermal transfer in optically
46
#' stimulated luminescence of quartz. Journal of Physics D: Applied Physics 40, 998-1006.
47
#'
48
#' Pagonis, V., Wintle, A.G., Chen, R., Wang, X.L., 2008. A theoretical model for a new dating protocol
49
#' for quartz based on thermally transferred OSL (TT-OSL).
50
#' Radiation Measurements 43, 704-708.
51
#'
52
#' @seealso \code{\link{model_LuminescenceSignals}}, \code{\link{translate_Sequence}}
53
#'
54
#' @examples
55
#'
56
#' #so far no example available
57
#'
58
#' @noRd
59
.simulate_RF <- function(
60
  temp,
61
  dose,
62
  dose_rate,
63
  RLumModel_ID = NULL,
64
  n,
65
  parms
66
){
67

68
# check input arguments ---------------------------------------------------
69

70
  ##check if temperature is > 0 K (-273 degree celsius)
71 6
  if(temp < -273){
72 6
    stop("\n [.simulate_RL()] Argument 'temp' has to be > 0 K!")
73
  }
74
  ##check if dose_rate is a positive number
75 6
  if(dose_rate < 0){
76 6
    stop("\n [.simulate_RL()] Argument 'dose_rate' has to be a positive number!")
77
  }
78

79
  ##check if dose is a positive number
80 6
  if(dose < 0){
81 6
    stop("\n [.simulate_RL()] Argument 'dose' has to be a positive number!")
82
  }
83

84
  ##check if n is a RLum object
85 6
  if(class(n) != "RLum.Results"){
86 6
    n <- n
87
  } else {
88 6
    n <- n$n
89
  }
90

91
# Set parameters for ODE ---------------------------------------------------
92

93
    ##============================================================================##
94
  # SETTING PARAMETERS FOR IRRADIATION
95
  #
96
  # R: electron-hole-production-rate (in Bailey 2004: 2.5e10, else: 5e7)
97
  # P: Photonflux (in Bailey 2004: wavelength [nm])
98
  # b: heating rate [deg. C/s]
99
  ##============================================================================##
100
  ## check if R is given in customized parameter sets
101 6
  if("R" %in% names(parms) && parms$R != 0){
102
    
103 6
    R <- dose_rate*parms$R
104
    
105
  } else {
106
    
107 6
    if(parms$model == "Bailey2004"){
108 6
      R <- dose_rate*2.5e10
109
    } else {
110
      
111 6
      if(parms$model == "Bailey2002"){
112 6
        R <- dose_rate*3e10
113
      } else {
114
        
115 6
        if(parms$model == "Friedrich2018"){
116 6
          R <- dose_rate*6.3e7
117
        } else {
118
          
119 6
          R <- dose_rate*5e7  # all other simulations
120
        }
121
      }
122
    }
123
  }
124 6
  P <- 0
125 6
  b <- 0
126

127
  ##============================================================================##
128
  # SETTING PARAMETERS FOR ODE
129
  ##============================================================================##
130

131 6
  times   <- seq(0, dose/(dose_rate), by = (dose/dose_rate)/1000)
132 6
  parameters.step <- .extract_pars(parameters.step = list(
133 6
    R = R,
134 6
    P = P,
135 6
    temp = temp,
136 6
    b = b,
137 6
    times = times,
138 6
    parms = parms))
139
  
140 6
  if(dose != 0){
141
  
142
  ##============================================================================##
143
  # SOLVING ODE (deSolve requiered)
144
  ##============================================================================##
145 6
  out <- deSolve::lsoda(y = n, times = times, parms = parameters.step, func = .set_ODE_Rcpp, rtol = 1e-4, atol = 1e-4,maxsteps = 10000)
146

147
  ##============================================================================##
148
  # CALCULATING RESULTS FROM ODE SOLVING
149
  ##============================================================================##
150

151 6
  signal <- .calc_signal(object = out, parameters = parameters.step)
152

153
  ##============================================================================##
154
  # CALCULATING CONCENTRATIONS FROM ODE SOLVING
155
  ##============================================================================##
156

157 6
  name <- c("RF")
158 6
  concentrations <- .calc_concentrations(
159 6
    data = out,
160 6
    times = times,
161 6
    name = name,
162 6
    RLumModel_ID = RLumModel_ID)
163

164
  ##============================================================================##
165
  # TAKING THE LAST LINE OF "OUT" TO COMMIT IT TO THE NEXT STEP
166
  ##============================================================================##
167

168 6
  return(Luminescence::set_RLum(class = "RLum.Results",
169 6
                                data = list(
170 6
                                n = out[length(times), -1],
171 6
                                RF.data = Luminescence::set_RLum(
172 6
                                  class = "RLum.Data.Curve",
173 6
                                  data = matrix(data = c(times, signal), ncol = 2),
174 6
                                  recordType = "RF",
175 6
                                  curveType = "simulated",
176 6
                                  info = list(
177 6
                                    curveDescripter = NA_character_),
178 6
                                    .pid = as.character(RLumModel_ID)
179
                                  ),
180 6
                                  temp = temp,
181 6
                                  concentrations = concentrations)
182
                  )
183
         )
184

185
  } else { ## dose == 0
186
    
187 6
    return(Luminescence::set_RLum(class = "RLum.Results",
188 6
                                  data = list(
189 6
                                    n = n,
190 6
                                    RF.data = Luminescence::set_RLum(
191 6
                                      class = "RLum.Data.Curve",
192 6
                                      data = matrix(data = c(times, 0), ncol = 2),
193 6
                                      recordType = "RF",
194 6
                                      curveType = "simulated",
195 6
                                      info = list(
196 6
                                        curveDescripter = NA_character_
197
                                        ),
198 6
                                      .pid = as.character(RLumModel_ID)
199
                                    ),
200 6
                                    temp = temp,
201 6
                                    concentrations = NULL)
202
                                  )
203
    )
204
    
205
    
206
  }
207
}

Read our documentation on viewing source code .

Loading