R-Lum / RLumModel

Compare 6909863 ... +10 ... bbe7e21

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.

Showing 5 of 33 files from the diff.

@@ -21,7 +21,7 @@
Loading
21 21
#'
22 22
#' @param model [character] (\bold{required}): set model to be used.
23 23
#' Available models are: `"Bailey2001"`, `"Bailey2002"`, `"Bailey2004"`, `"Pagonis2007"`,
24 -
#' `"Pagonis2008"`, `"Friedrich2017"`, `"Friedrich2018"`. If model is indeed missing,
24 +
#' `"Pagonis2008"`, `"Friedrich2017"`, `"Friedrich2018"`, `"Peng2022"`. If model is indeed missing,
25 25
#' the list of allowed keywords is returned.
26 26
#'
27 27
#' @return This function returns a [list] with all necessary parameters for
@@ -62,6 +62,8 @@
Loading
62 62
#' for quartz based on thermally transferred OSL (TT-OSL).
63 63
#' Radiation Measurements 43, 704-708.
64 64
#'
65 +
#' Peng, J., Wang, X., Adamiec, G., 2022. The build-up of the laboratory-generated dose-response curve and underestimation of equivalent dose for quartz OSL in the high dose region: A critical modelling study. Quaternary Geochronology 67, 101231.
66 +
#'
65 67
#' @examples
66 68
#'
67 69
#' pars <- .set_pars("Bailey2001")
@@ -80,6 +82,7 @@
Loading
80 82
    "Bailey2002",
81 83
    "Friedrich2017",
82 84
    "Friedrich2018",
85 +
    "Peng2022",
83 86
    "customized",
84 87
    "customised"
85 88
  )
@@ -283,6 +286,26 @@
Loading
283 286
      model = model
284 287
    ),
285 288
289 +
    Peng2022 = list(
290 +
      N = c(1.5e7, 1e7, 1e9, 2.5e8, 5e10, 1.65e8, 5e9, 5e9, 1e11),
291 +
      E = c(0.97, 1.55, 1.7, 1.72, 2, 1.41, 1.65, 5, 5),
292 +
      s = c(5e12, 5e14, 1e13, 1e14, 1e10, 5e13, 5e14, 1e13, 1e13),
293 +
      A = c(1e-8, 1e-8, 1e-9, 5e-10, 3.33e-11, 5e-7, 1e-9, 1e-10, 1e-9),
294 +
      B = c(0, 0, 0, 0, 0, 5e-9, 5e-10, 1e-10, 1e-10),
295 +
      Th = c(0.75, 0, 6, 4.5, 0),
296 +
      E_th = c(0.1, 0, 0.1, 0.13, 0),
297 +
      n =  set_RLum(class = "RLum.Results", data = list(
298 +
        n = c(9.169767e-03, 7.619894e+04, 1.291564e+08, 7.432290e+06, 2.690423e+10,
299 +
              5.741230e+06, 6.779304e+07, 1.591234e+08, 2.680824e+10, 2.450977e-07, 4.263486e-07),
300 +
        temp = 20)),
301 +
      k_B = k_B,
302 +
      W = W,
303 +
      K = K,
304 +
      units  = units,
305 +
      names = names,
306 +
      model = model
307 +
    ),
308 +
286 309
    customized = list(
287 310
      n =  set_RLum(class = "RLum.Results", data = list(
288 311
        n = rep(0,4), temp = 20, model = model)),
@@ -306,7 +329,6 @@
Loading
306 329
    )
307 330
  )
308 331
309 -
310 332
  switch(model,
311 333
        "Bailey2001" = {
312 334
          return(parameter.list$Bailey2001)
@@ -336,6 +358,10 @@
Loading
336 358
          return(parameter.list$Friedrich2018)
337 359
        },
338 360
361 +
        "Peng2022" = {
362 +
          return(parameter.list$Peng2022)
363 +
        },
364 +
339 365
        "customized" = {
340 366
          return(parameter.list$customized)
341 367
        },

@@ -70,7 +70,7 @@
Loading
70 70
#' Risø sequence editor. To simulate SAR measurements there is an extra option to set the sequence list (cf. details).
71 71
#
72 72
#' @param model [character] (**required**): set model to be used. Available models are:
73 -
#' `"Bailey2001"`, `"Bailey2002"`, `"Bailey2004"`, `"Pagonis2007"`, `"Pagonis2008"`, `"Friedrich2017"`, `"Friedrich2018"` and for own models `"customized"` (or `"customised"`).
73 +
#' `"Bailey2001"`, `"Bailey2002"`, `"Bailey2004"`, `"Pagonis2007"`, `"Pagonis2008"`, `"Friedrich2017"`, `"Friedrich2018"`, `"Peng2022`" and for own models `"customized"` (or `"customised"`).
74 74
#' Note: When model = `"customized"`/`"customised` is set, the argument `own_parameters` has to be set.
75 75
#'
76 76
#' @param lab.dose_rate [numeric] (*with default*): laboratory dose rate in XXX
@@ -136,6 +136,10 @@
Loading
136 136
#' Pagonis, V., Lawless, J., Chen, R., Anderson, C., 2009. Radioluminescence in Al2O3:C - analytical and numerical
137 137
#' simulation results. Journal of Physics D: Applied Physics 42, 175107 (9pp).
138 138
#'
139 +
#' Peng, J., Wang, X., Adamiec, G., 2022. The build-up of the laboratory-generated dose-response curve and
140 +
#' underestimation of equivalent dose for quartz OSL in the high dose region: A critical modelling study.
141 +
#' Quaternary Geochronology 67, 101231.
142 +
#'
139 143
#' Soetaert, K., Cash, J., Mazzia, F., 2012. Solving differential equations in R.
140 144
#' Springer Science & Business Media.
141 145
#'

@@ -0,0 +1,203 @@
Loading
1 +
#'@title Trace parameter state evolution
2 +
#'
3 +
#'@description Traces the evolution of the concentrations in the different
4 +
#'levels over different simulation steps. For instance, a sequence consisting
5 +
#'of one TL and one OSL step has two iterations. For each step the end concentration
6 +
#'is extracted. This way, the evolution of the system can be traced throughout
7 +
#'a sequence.
8 +
#'
9 +
#'@param object [Luminescence::RLum.Analysis-class] (**required**): input object
10 +
#'created by the function [model_LuminescenceSignals]. The input can be a list of such
11 +
#'objects
12 +
#'
13 +
#'@param step [character] (*optional*): filter the input object to pick particular
14 +
#'steps, the input is passed to [Luminescence::get_RLum]
15 +
#'
16 +
#'@param plot [logical] (*with default*): enables/disables plot output
17 +
#'
18 +
#'@param ... optional arguments to be passed to control the plot output. Supported
19 +
#'are `xlab`, `ylab`, `col`, `type`, `bg`, `main`. Where meaningful, parameters can be provided
20 +
#'as vectors. Vectors short than the number of plots are recycled.
21 +
#'
22 +
#'@return Returns a plot and [list] with [matrix] objects of the parameter evolution. If
23 +
#'`object` is a [list] the output is a nested [list]
24 +
#'
25 +
#'@section Function version: 0.1.0
26 +
#'
27 +
#'@author Sebastian Kreutzer, Geography & Earth Sciences, Aberystwyth University (United Kingdom)
28 +
#'
29 +
#'@examples
30 +
#'
31 +
#' sequence <-
32 +
#' list(
33 +
#' IRR = c(20, 10, 1),
34 +
#' TL = c(20, 400, 5),
35 +
#' IRR = c(20, 10, 1),
36 +
#' TL = c(20, 400, 5))
37 +
#'
38 +
#'##model sequence
39 +
#'model.output <- model_LuminescenceSignals(
40 +
#'  sequence = sequence,
41 +
#'  verbose = FALSE,
42 +
#'  plot = FALSE,
43 +
#'  model = "Bailey2001")
44 +
#'
45 +
#'## trace
46 +
#'trace_ParameterStateEvolution(model.output)
47 +
#'
48 +
#'@md
49 +
#'@export
50 +
trace_ParameterStateEvolution <- function(
51 +
  object,
52 +
  step = NULL,
53 +
  plot = TRUE,
54 +
  ...
55 +
){
56 +
57 +
# Self call ---------------------------------------------------------------
58 +
  if (is(object, "list")) {
59 +
    arguments <- list(step = step[1], plot = plot[1], list(...))
60 +
61 +
    return(lapply(object, function(x){
62 +
      do.call("trace_ParameterStateEvolution", c(list(object = x), arguments))
63 +
    }))
64 +
65 +
  }
66 +
67 +
# Check incoming ----------------------------------------------------------
68 +
  if (!is(object, "RLum.Analysis"))
69 +
    stop("[trace_ParameterStateEvolution()] object is not of class 'RLum.Analysis!",
70 +
         call. = FALSE)
71 +
72 +
  if (object@originator != "model_LuminescenceSignals")
73 +
    stop("[trace_ParameterStateEvolution()] object was not produced by model_LuminescencerSignals()!",
74 +
         call. = FALSE)
75 +
76 +
  if(!any(grepl("conc.", names(object))))
77 +
    stop("[trace_ParameterStateEvolution()] No concentration record found, did you subset your object already?",
78 +
          call. = FALSE)
79 +
80 +
# Prepare data ------------------------------------------------------------
81 +
  if (!is.null(step))
82 +
    object <- get_RLum(object, recordType = step[1], drop = FALSE)
83 +
84 +
  if (length(object) == 0)
85 +
    stop("[trace_ParameterStateEvolution()] object has length zero!", call. = FALSE)
86 +
87 +
# Extract parameter evolution ---------------------------------------------
88 +
  ## get list of parameters ... but only the names of the concentration
89 +
  param_names <- names(object)
90 +
  param_names <- unique(unlist(
91 +
    regmatches(
92 +
      x = param_names,
93 +
      m = regexec("conc\\.\\s.+(?=\\s\\()", param_names, perl = TRUE))))
94 +
95 +
  ## extract start and end values of each state
96 +
  m_list <- lapply(param_names, function(p){
97 +
    vapply(get_RLum(object, recordType = p), function(x){
98 +
          x@data[c(1,nrow(x@data)),2]
99 +
        }, numeric(2))
100 +
101 +
  })
102 +
103 +
  ## assign rownames
104 +
  m_list <- lapply(m_list, function(x){
105 +
    rownames(x) <- c("start", "end")
106 +
    x
107 +
  })
108 +
109 +
  ## assigning level names
110 +
  names(m_list) <- param_names
111 +
112 +
# Plotting ----------------------------------------------------------------
113 +
  if (plot) {
114 +
    ## obtain number of plots to produce
115 +
    n <- length(m_list)
116 +
117 +
    ## plot area
118 +
    if (n > 1) {
119 +
      op <- par(
120 +
        mfrow = c(n, 1),
121 +
        mar = c(0, 5, 0, 5),
122 +
        omi = c(0.5, 0, 0.5, 0)
123 +
      )
124 +
      on.exit(par(op))
125 +
    }
126 +
127 +
    ## set plot settings
128 +
    plot_settings <- modifyList(list(
129 +
      main = "State parameter evolution (end concentrations)",
130 +
      xlab = "Evolution index",
131 +
      ylab = regmatches(param_names, regexec("[^conc\\.].+", param_names)),
132 +
      col = khroma::colour("muted")(),
133 +
      bg = TRUE,
134 +
      type = "l"
135 +
136 +
    ), list(...))
137 +
138 +
    ## expand parameters, this saves us a lot of trouble later on
139 +
    plot_settings <- lapply(plot_settings, rep, length.out = n)
140 +
141 +
    ## print single plots
142 +
    for (i in 1:n) {
143 +
      ## odd and even change of plots
144 +
      if (i %% 2 != 0) {
145 +
        plot(NA, NA,
146 +
          xaxt = "n",
147 +
          xlab = if (n > 1) "" else plot_settings$xlab[1],
148 +
          xlim = c(1,ncol(m_list[[i]])),
149 +
          ylim = range(m_list[[i]]),
150 +
          ylab = plot_settings$ylab[i],
151 +
          frame = FALSE)
152 +
153 +
        ## add background
154 +
        if (plot_settings$bg[1]) {
155 +
          rect(
156 +
            par("usr")[1],
157 +
            par("usr")[3],
158 +
            par("usr")[2],
159 +
            par("usr")[4],
160 +
            col = rgb(0, 0, 0, 0.1),
161 +
            border = FALSE
162 +
          )
163 +
        }
164 +
165 +
        ## add lines
166 +
        lines(
167 +
          x = 1:ncol(m_list[[i]]),
168 +
          y = m_list[[i]][2,],
169 +
          type = plot_settings$type[i],
170 +
          col = plot_settings$col[i])
171 +
172 +
        ## add title
173 +
        if (i == 1)
174 +
          mtext(side = 3, plot_settings$main[1], outer = TRUE, line = 1)
175 +
176 +
      } else {
177 +
        plot(
178 +
          x = 1:ncol(m_list[[i]]),
179 +
          y = m_list[[i]][2,],
180 +
          xaxt = "n",
181 +
          yaxt = "n",
182 +
          type = plot_settings$type[i],
183 +
          frame = FALSE,
184 +
          ylab = "",
185 +
          col = plot_settings$col[i])
186 +
        axis(side = 4)
187 +
        mtext(side = 4, text = plot_settings$ylab[i], line = 3, cex = 0.7)
188 +
      }
189 +
190 +
    }
191 +
192 +
    ## add x-axis
193 +
    axis(side = 1, at = pretty(1:ncol(m_list[[1]])), las = 1)
194 +
195 +
    if (n > 1)
196 +
      mtext(side = 1, text = plot_settings$xlab[1], cex = 0.7, line = 2)
197 +
  }
198 +
199 +
# Return ------------------------------------------------------------------
200 +
 return(m_list)
201 +
202 +
}
203 +

@@ -1,6 +1,6 @@
Loading
1 -
#' sequence step irradiation
1 +
#' @title sequence step irradiation
2 2
#'
3 -
#' This function simulates the irradiaton of quartz in the energy-band-model.
3 +
#' @description This function simulates the irradiation of quartz in the energy-band-model.
4 4
#'
5 5
#' @param temp \code{\link{cnumeric}} (\bold{required}): temperature [deg. C] at which the dose should be applied
6 6
#'
@@ -58,16 +58,15 @@
Loading
58 58
  }
59 59
60 60
  ##check if n is a RLum object
61 -
  if(class(n) != "RLum.Results"){
61 +
  if(class(n)[1] != "RLum.Results"){
62 62
    n <- n
63 63
  } else {
64 64
    n <- n$n
65 65
  }
66 66
67 67
# solving ODE ---------------------------------------------------
68 -
  
69 68
  if(dose != 0){
70 -
  
69 +
71 70
  ##============================================================================##
72 71
  # SETTING PARAMETERS FOR IRRADIATION
73 72
  #
@@ -77,27 +76,19 @@
Loading
77 76
  ##============================================================================##
78 77
  ## check if R is given in customized parameter sets
79 78
  if("R" %in% names(parms) && parms$R != 0){
80 -
    
81 79
    R <- dose_rate*parms$R
82 -
    
80 +
83 81
  } else {
82 +
    R <- dose_rate * 5e7  # all other simulations
83 +
84 +
    if(parms$model == "Bailey2004")
85 +
      R <- dose_rate * 2.5e10
84 86
85 -
    if(parms$model == "Bailey2004"){
86 -
      R <- dose_rate*2.5e10
87 -
    } else {
88 -
      
89 -
      if(parms$model == "Bailey2002"){
90 -
        R <- dose_rate*3e10
91 -
      } else {
92 -
        
93 -
        if(parms$model == "Friedrich2018"){
94 -
          R <- dose_rate*6.3e7
95 -
        } else {
96 -
        
97 -
          R <- dose_rate*5e7  # all other simulations
98 -
        }
99 -
      }
100 -
    }
87 +
    if(parms$model == "Bailey2002")
88 +
      R <- dose_rate * 3e10
89 +
90 +
    if(parms$model == "Friedrich2018")
91 +
      R <- dose_rate * 6.3e7
101 92
  }
102 93
103 94
  P <- 0
@@ -106,7 +97,6 @@
Loading
106 97
  ##============================================================================##
107 98
  # SETTING PARAMETERS FOR ODE
108 99
  ##============================================================================##
109 -
110 100
  times   <- seq(0, dose/(dose_rate), by = (dose/dose_rate)/100)
111 101
  parameters.step <- .extract_pars(parameters.step = list(
112 102
    R = R,
@@ -115,31 +105,32 @@
Loading
115 105
    b = b,
116 106
    times = times,
117 107
    parms = parms))
118 -
  
108 +
119 109
  ##============================================================================##
120 -
  # SOLVING ODE (deSolve requiered)
110 +
  # SOLVING ODE (deSolve required)
121 111
  ##============================================================================##
122 -
  
123 -
  out <- deSolve::lsoda(y = n, times = times, parms = parameters.step, func = .set_ODE_Rcpp, rtol = 1e-6, atol = 1e-6);
112 +
  out <- deSolve::lsoda(
113 +
    y = n,
114 +
    times = times,
115 +
    parms = parameters.step,
116 +
    func = .set_ODE_Rcpp,
117 +
    rtol = 1e-6,
118 +
    atol = 1e-6);
124 119
125 120
  ##============================================================================##
126 121
  # TAKING THE LAST LINE OF "OUT" TO COMMIT IT TO THE NEXT STEP
127 122
  ##============================================================================##
128 -
129 -
  return(Luminescence::set_RLum(class = "RLum.Results",
130 -
                  data = list(
131 -
                    n = out[length(times),-1],
132 -
                    temp = temp
133 -
                  )))
123 +
  return(Luminescence::set_RLum(
124 +
    class = "RLum.Results",
125 +
    data = list(
126 +
    n = out[length(times),-1],temp = temp)))
134 127
135 128
  } else {
136 -
    
137 -
  return(Luminescence::set_RLum(class = "RLum.Results",
138 -
                                data = list(
139 -
                                  n = n,
140 -
                                  temp = temp
141 -
                                )))
142 -
  
143 -
  
144 -
}
129 +
  return(Luminescence::set_RLum(
130 +
    class = "RLum.Results",
131 +
    data = list(
132 +
    n = n,
133 +
    temp = temp)))
134 +
135 +
 }
145 136
}

@@ -388,7 +388,7 @@
Loading
388 388
  class = "RLum.Analysis",
389 389
  records = output.model,
390 390
  protocol = model,
391 -
  originator = "model_LuminescenceSignals()",
391 +
  originator = "model_LuminescenceSignals",
392 392
  info = list(
393 393
    sequence = sequence,
394 394
    parms = parms,

Learn more Showing 1 files with coverage changes found.

New file R/trace_ParameterStateEvolution.R
New
Loading file...
Files Coverage
R 0.74% 92.66%
src 100.00%
Project Totals (21 files) 93.01%
Loading