R-Lum / RLumModel
1
#' sequence step heating/cooling between different simulation steps
2
#'
3
#' This function simulates the heating/cooling 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_begin \code{\link{numeric}} (\bold{required}): endtemperature [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 n \code{\link{numeric}} or \code{\linkS4class{RLum.Results}} (\bold{required}):
12
#' concentration of electron-/holetraps, valence- and conduction band
13
#' from step before. This is necessary to get the boundary condition for the ODEs.
14
#'
15
#' @param parms \code{\linkS4class{RLum.Results}} (\bold{required}): The specific model parameters are used to simulate
16
#' numerical quartz luminescence results.
17
#'
18
#' @param \dots further arguments and graphical parameters passed to
19
#' \code{\link{plot.default}}. See details for further information
20
#'
21
#' @return This function returns an RLum.Results object from the heating/cooling simulation.
22
#'
23
#' @section Function version: 0.1.1
24
#'
25
#' @author Johannes Friedrich, University of Bayreuth (Germany),
26
#'
27
#' @references
28
#'
29
#' Bailey, R.M., 2001. Towards a general kinetic model for optically and thermally stimulated
30
#' luminescence of quartz. Radiation Measurements 33, 17-45.
31
#'
32
#' @seealso \code{\link{simulate_TL}}
33
#'
34
#' @examples
35
#'
36
#' #so far no example available
37
#'
38
#' @noRd
39
.simulate_heating <- function(
40
  temp_begin,
41
  temp_end,
42
  heating_rate,
43
  n,
44
  parms
45
){
46

47
# check input arguments ---------------------------------------------------
48

49
  ##check if heatingrate has the rigth algebraic sign
50 5
  if((temp_begin < temp_end && heating_rate < 0)||(temp_begin > temp_end & heating_rate > 0)){
51 5
    stop("\n [.simulate_heating()] Heatingrate has the wrong algebraic sign!")
52
  }
53

54
  ##check if temperature is > 0 K (-273 degree celsius)
55 5
  if(temp_begin < -273 ||temp_end < -273){
56 5
    stop("\n [.simulate_heating()] Argument 'temp' has to be > 0 K!")
57
  }
58

59
  ##check if object is of class RLum.Results
60 5
  if(class(n) != "RLum.Results"){
61 5
    n <- n
62
  } else {
63 5
    n <- n$n
64
  }
65

66
# Set parameters for ODE ---------------------------------------------------
67

68
  ##============================================================================##
69
  # SETTING PARAMETERS FOR HEATING
70
  #
71
  # R: electron-hole-production-rate (in Bailey 2004: 2.5e10, else: 5e7) = 0
72
  # P: Photonflux (in Bailey 2004: wavelength [nm]) = 0
73
  # b: heating rate [deg. C/s]
74
  ##============================================================================##
75

76 5
  R <- 0
77 5
  P <- 0
78 5
  b <- heating_rate
79

80
  ##============================================================================##
81
  # SETTING PARAMETERS FOR ODE
82
  ##============================================================================##
83

84 5
  times <- seq(from = 0, to = (temp_end-temp_begin)/b, by = 0.1)
85 5
  parameters.step <- .extract_pars(parameters.step = list(
86 5
    R = R,
87 5
    P = P,
88 5
    temp = temp_begin,
89 5
    b = b,
90 5
    times = times,
91 5
    parms = parms))
92
  ##============================================================================##
93
  # SOLVING ODE (deSolve requiered)
94
  ##============================================================================##
95 5
  out <- deSolve::lsoda(y = n, times = times, parms = parameters.step, func = .set_ODE_Rcpp, maxsteps = 100000, rtol = 1e-4, atol = 1e-4)
96
  ##============================================================================##
97

98
  ##============================================================================##
99
  # TAKING THE LAST LINE OF "OUT" TO COMMIT IT TO THE NEXT STEP
100
  ##============================================================================##
101

102 5
  return(Luminescence::set_RLum(class = "RLum.Results",
103 5
                  data = list(
104 5
                    n = out[length(times),-1],
105 5
                    temp = temp_end
106
                  )))
107

108
}#end function

Read our documentation on viewing source code .

Loading