R-Lum / RLumModel
1
#' sequence step heating/cooling between different simulation steps
2
#'
3
#' This function simulates the heating/cooling and simultaneously irradiation and signal detection  of quartz in the energy-band-model.
4
#'
5
#' @param temp_begin \code{\link{numeric}} (\bold{required}): initial temperature [deg. C] of the TL-simulation
6
#'
7
#' @param temp_end \code{\link{numeric}} (\bold{required}): end temperature [deg. C] of the TL-simulation
8
#'
9
#' @param heating_rate \code{\link{numeric}} (\bold{required}): heatingrate in [deg. C/s] or [K/s]
10
#' 
11
#' @param dose_rate \code{\link{numeric}} (\bold{required}): dose-rate in [Gy/s] for the irradiation during heating
12
#'
13
#' @param RLumModel_ID \code{\link{numeric}} (optional): A ID-number for the TL-step. This ID
14
#' is pass down to \link{calc_concentrations} so all concentrations had the same ID as the
15
#' sequence step they were calculated from. This ID is identic to the sequence step in "sequence".
16
#' 
17
#' @param n \code{\link{numeric}} or \code{\linkS4class{RLum.Results}} (\bold{required}):
18
#' concentration of electron-/holetraps, valence- and conduction band
19
#' from step before. This is necessary to get the boundary condition for the ODEs.
20
#'
21
#' @param parms \code{\linkS4class{RLum.Results}} (\bold{required}): The specific model parameters are used to simulate
22
#' numerical quartz luminescence results.
23
#'
24
#' @return This function returns an RLum.Results object from the heating/cooling simulation during irradiation.
25
#'
26
#' @section Function version: 0.1.0 [2017-09-01]
27
#'
28
#' @author Johannes Friedrich, University of Bayreuth (Germany)
29
#'
30
#' @references
31
#'
32
#' Bailey, R.M., 2001. Towards a general kinetic model for optically and thermally stimulated
33
#' luminescence of quartz. Radiation Measurements 33, 17-45.
34
#' 
35
#' Poolton, N. R. J., et al., 2001. An automated system for the analysis of variable temperature radioluminescence.
36
#' Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 179.4, 575-584.
37
#' 
38
#' Wintle, A. G. 1975. Thermal quenching of thermoluminescence in quartz. Geophysical Journal International, 41(1), 107-113.#' 
39
#' 
40
#'
41
#' @seealso \code{\link{simulate_TL}}, code{\link{simulate_RF}}
42
#'
43
#' @examples
44
#'
45
#' #so far no example available
46
#'
47
#' @noRd
48
.simulate_RF_and_heating <- function(
49
  temp_begin,
50
  temp_end,
51
  heating_rate,
52
  dose_rate,
53
  RLumModel_ID = NULL,
54
  n,
55
  parms
56
){
57
  
58
  # check input arguments ---------------------------------------------------
59
  
60
  ##check if heatingrate has the rigth algebraic sign
61 6
  if((temp_begin < temp_end && heating_rate < 0)||(temp_begin > temp_end & heating_rate > 0)){
62 6
    stop("\n [.simulate_RF_and_heating()] Heatingrate has the wrong algebraic sign!")
63
  }
64
  
65
  ##check if temperature is > 0 K (-273 degree celsius)
66 6
  if(temp_begin < -273 ||temp_end < -273){
67 6
    stop("\n [.simulate_RF_and_heating()] Argument 'temp' has to be > 0 K!")
68
  }
69
  
70
  ##check if dose_rate is a positive number
71 6
  if(dose_rate < 0){
72 6
    stop("\n [.simulate_RF_and_heating()] Argument 'dose_rate' has to be a positive number!")
73
  }
74
  
75
  ##check if object is of class RLum.Results
76 6
  if(class(n) != "RLum.Results"){
77 6
    n <- n
78
  } else {
79 6
    n <- n$n
80
  }
81
  
82
  # Set parameters for ODE ---------------------------------------------------
83
  
84
    
85
    ##============================================================================##
86
    # SETTING PARAMETERS FOR IRRADIATION
87
    #
88
    # R: electron-hole-production-rate (in Bailey 2004: 2.5e10, Bailey 2002: 3e10, else: 5e7)
89
    # P: Photonflux (in Bailey 2004: wavelength [nm])
90
    # b: heating rate [deg. C/s]
91
    ##============================================================================##
92
    ## check if R is given in customized parameter sets
93 6
    if("R" %in% names(parms) && parms$R != 0){
94
      
95 0
      R <- dose_rate*parms$R
96
      
97
    } else {
98
      
99 6
      if(parms$model == "Bailey2004"){
100 6
        R <- dose_rate*2.5e10
101
      } else {
102
        
103 6
        if(parms$model == "Bailey2002"){
104 6
          R <- dose_rate*3e10
105
        } else {
106
          
107 6
          if(parms$model == "Friedrich2018"){
108 6
            R <- dose_rate*6.3e7
109
          } else {
110
            
111 6
            R <- dose_rate*5e7  # all other simulations
112
          }
113
        }
114
      }
115
    }
116
    
117 6
    P <- 0
118 6
    b <- heating_rate
119
    
120
    ##============================================================================##
121
    # SETTING PARAMETERS FOR ODE
122
    ##============================================================================##
123
    
124 6
    times <- seq(from = 0, to = (temp_end-temp_begin)/b, by = 0.1)
125 6
    parameters.step <- .extract_pars(parameters.step = list(
126 6
      R = R,
127 6
      P = P,
128 6
      temp = temp_begin,
129 6
      b = b,
130 6
      times = times,
131 6
      parms = parms))
132
    
133
    ##============================================================================##
134
    # SOLVING ODE (deSolve requiered)
135
    ##============================================================================##
136
    
137 6
    out <- deSolve::lsoda(y = n, times = times, parms = parameters.step, func = .set_ODE_Rcpp, rtol = 1e-6, atol = 1e-6);
138
    
139
    ##============================================================================##
140
    # CALCULATING RESULTS FROM ODE SOLVING
141
    ##============================================================================##
142
    
143
    
144 6
    signal <- .calc_signal(object = out, parameters = parameters.step)
145 6
    TSkala <- times*b+temp_begin
146
    
147
    ##============================================================================##
148
    # CALCULATING CONCENTRATIONS FROM ODE SOLVING
149
    ##============================================================================##
150
    
151 6
    name <- c("RF_heating")
152 6
    concentrations <- .calc_concentrations(
153 6
      data = out,
154 6
      times = TSkala,
155 6
      name = name,
156 6
      RLumModel_ID = RLumModel_ID)
157
    
158
    ##============================================================================##
159
    # TAKING THE LAST LINE OF "OUT" TO COMMIT IT TO THE NEXT STEP
160
    ##============================================================================##
161
    
162 6
    return(Luminescence::set_RLum(class = "RLum.Results",
163 6
                                  data = list(
164 6
                                    n = out[length(times),-1] ,
165 6
                                    RF_heating.data = Luminescence::set_RLum(
166 6
                                      class = "RLum.Data.Curve",
167 6
                                      data = matrix(data = c(TSkala, signal),ncol = 2),
168 6
                                      recordType = "RF_heating",
169 6
                                      curveType = "simulated",
170 6
                                      info = list(
171 6
                                        curveDescripter = NA_character_),
172 6
                                      .pid = as.character(RLumModel_ID)
173
                                    ),
174 6
                                    temp = temp_end,
175 6
                                    concentrations = concentrations)
176
    )
177
    )
178
    
179
}#end function

Read our documentation on viewing source code .

Loading