R/use_DRAC4flint.R

In TLdating: Tools for Thermoluminescences Dating

#' Calculation of the dose rate for a flint
#'
#' This function allows to estimating the dose rate for a burnt flint object which was surrounded by sediment.
#' It call the \link{use_DRAC} and \link{calc_CosmicDoseRate} functions from the R package '\link{Luminescence}'.
#' The function 'use_DRAC' is only compatible wth DRAC version 1.1.
#'
#' @param data
#'  \link{list}: data object create throught the function \link{template_DRAC4flint}.
#'
#' @param notification
#'  \link{logical} (default): set to \code{FALSE} if you don't want to see the legal notification.
#'
#' @return
#'  This function return a \linkS4class{TLum.Results} object containing the Age estimation, the dose rates
#'  (total, internal, external, environmental, alpha, beta, gamma and cosmic), the equivalent dose used and their uncertainties.
#'
#' @author David Strebler
#'
#' @export use_DRAC4flint

use_DRAC4flint <- function (

  data,

  notification = TRUE
){
  ## ------------------------------------------------------------------------------------------

  if(!is(data,"DRAC4flint.list")){
    stop("[use_DRAC4flint] Error: data has to be create using the 'template_DRAC4flint' function.")
  }

  ## ------------------------------------------------------------------------------------------

  ## LEGAL NOTICE ----

  legal.notice <- paste("\t	-------------------- IMPORTANT NOTE ------------------------ \n",
                        "\t This function relies on the functions 'use_DRAC' and 'calc_CosmicDoseRate' \n",
                        "\t from the R package 'Luminescence' to estimate the different dose rate. \n",
                        "\t The function 'use_DRAC' is only compatible with DRAC version 1.1. \n",
                        "\t Before using this function make sure that this is the correct version, otherwise expect unspecified errors.\n",
                        "\t Please ensure you cite the use of DRAC in your work, published or otherwise. Please cite the website name and  \n",
                        "\t version (e.g. DRAC v1.1) and the accompanying journal article:  \n",
                        "\t Durcan, J.A., King, G.E., Duller, G.A.T., 2015. DRAC: Dose rate and age calculation for trapped charge  \n",
                        "\t dating. Quaternary Geochronology 28, 54-61. \n",
                        "\t Please ensure you also cite the use of the R package 'Luminescence' (cf. citation('Luminescence')). \n")

  if(notification){
    message(legal.notice)
  }

  comment <- list(legal=legal.notice)

  project <- data$info$project
  sample <- data$info$sample
  mineral <- data$info$mineral
  conversion.factors <- data$info$conversion.factors

  int.alpha.attenuation <- data$info$alpha.size.attenuation
  int.beta.attenuation <- data$info$beta.size.attenuation
  int.etch.beta.attenuation <- data$info$beta.etch.attenuation

  date <- data$info$date

  ## ------------------------------------------------------------------------------------------
  # Check values
  if(is.null(project)){
    project <- "Unknown"
    warning("[use_DRAC4flint] Warning: no project name")

  }else if(!is.character(project)){
    stop("[use_DRAC4flint] Error: data$info$project has to be of type 'character'.")
  }

  if(is.null(sample)){
    sample <- "Unknown"
    warning("[use_DRAC4flint] Warning: no sample name")

  }else if(!is.character(sample)){
    stop("[use_DRAC4flint] Error: data$info$sample has to be of type 'character'.")
  }

  if(is.null(mineral)){
    stop("[use_DRAC4flint] Error: mineral unknown")
  }else if(!is.character(mineral)){
    stop("[use_DRAC4flint] Error: data$info$mineral has to be of type 'character'.")
  }else if(!(mineral %in% c("Q", "F", "PM"))){
    stop("[use_DRAC4flint] Error: data$info$mineral can only be 'Q', 'F' or 'PM'.")
  }

  if(is.null(conversion.factors)){
    stop("[use_DRAC4flint] Error: no conversion factors selected")
  }else if(!is.character(conversion.factors)){
    stop("[use_DRAC4flint] Error: data$info$conversion.factors has to be of type 'character'.")
  }else if(!(conversion.factors %in% c("AdamiecAitken1998", "Guerinetal2011", "Liritzisetal2013"))){
    stop("[use_DRAC4flint] Error: data$info$conversion.factors can only be 'AdamiecAitken1998', 'Guerinetal2011' or 'Liritzisetal2013'.")
  }

  if(is.null(int.alpha.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$alpha.size.attenuation' is null.")
  }else if(!is.character(int.alpha.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$alpha.size.attenuation' is not of type 'character'.")
  }else if(!(int.alpha.attenuation %in% c("Bell1980","Brennanetal1991"))){
    stop("[use_DRAC4flint] Error: 'data$info$alpha.size.attenuation' can only be 'Bell1980' or 'Brennanetal1991'.")
  }

  if(is.null(int.beta.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$beta.size.attenuation' is null.")
  }else if(!is.character(int.beta.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$beta.size.attenuation' is not of type 'character'.")
  }else if(!(int.beta.attenuation %in% c("Mejdahl1979","Brennan2003", "Guerinetal2012-Q", "Guerinetal2012-F"))){
    stop("[use_DRAC4flint] Error: 'data$info$beta.size.attenuation' can only be 'Mejdahl1979', 'Brennan2003', 'Guerinetal2012-Q' or 'Guerinetal2012-F'.")
  }else if(int.beta.attenuation == "Guerinetal2012-Q" && mineral != "Q"){
    warning("[use_DRAC4flint]  Warning you use 'Guerinetal2012-Q' for 'data$info$beta.size.attenuation' but material is not 'Q'.")
  }else if(int.beta.attenuation == "Guerinetal2012-F" && mineral != "F"){
    warning("[use_DRAC4flint]  Warning you use 'Guerinetal2012-F' for 'data$info$beta.size.attenuation' but material is not 'F'.")
  }

  if(is.null(int.etch.beta.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$beta.etch.attenuation' is null.")
  }else if(!is.character(int.etch.beta.attenuation)){
    stop("[use_DRAC4flint] Error: 'data$info$beta.etch.attenuation' is not of type 'character'.")
  }else if(!(int.etch.beta.attenuation %in% c("Bell1979","Brennan2003"))){
    stop("[use_DRAC4flint] Error: 'data$info$beta.etch.attenuation' can only be 'Bell1979' or 'Brennan2003'.")
  }

  if(is.null(date)){
    date <- as.numeric(format(Sys.Date(),"%Y"))
  }else if(!is.numeric(date)){
    stop("[use_DRAC4flint] Error: 'data$info$date' is not of type 'numeric'.")
  }

  ## ------------------------------------------------------------------------------------------

  # Equivalent dose
  De <- data$De$De
  De.err <- data$De$De.err

  ## ------------------------------------------------------------------------------------------
  # Check values
  if(is.null(De) || is.null(De.err)){
    De <- 0
    De.err <- 0
    warning("[use_DRAC4flint] Warning: data$De$De or data$De$De.err is missing. De and De.err are considered as null")

  }else if(!is.numeric(De) || !is.numeric(De)){
    stop("[use_DRAC4flint] Error: data$De$De has to be of type 'numeric'.")
  }else if(De < 0 || De.err < 0){
    stop("[use_DRAC4flint] Error: data$De$De or data$De$De.err < 0.")
  }
  ## ------------------------------------------------------------------------------------------

  # Internal dose rate
  # By concentration
  int.U <- data$flint$Dr$U
  int.U.err <- data$flint$Dr$U.err
  int.Th <- data$flint$Dr$Th
  int.Th.err <- data$flint$Dr$Th.err
  int.K <- data$flint$Dr$K
  int.K.err <- data$flint$Dr$K.err
  int.Rb <- data$flint$Dr$Rb
  int.Rb.err <- data$flint$Dr$Rb.err

  int.K2Rb <- data$flint$Dr$K2Rb


  # Direct evaluation
  int.alpha <- data$flint$Dr$alpha
  int.alpha.err <- data$flint$Dr$alpha.err
  int.beta <- data$flint$Dr$beta
  int.beta.err <- data$flint$Dr$beta.err
  int.gamma <- data$flint$Dr$gamma
  int.gamma.err <- data$flint$Dr$gamma.err


  ## ------------------------------------------------------------------------------------------
  # Check values

  # Dose rate direct estimation
  if(is.null(int.alpha) || is.null(int.alpha.err)){
    int.alpha <- "X"
    int.alpha.err <- "X"
  }else if(!is.numeric(int.alpha) || !is.numeric(int.alpha.err)){
    stop("[use_DRAC4flint] Error: data$flint$Dr$alpha or data$flint$Dr$alpha.err is not of type 'numeric'.")
  }else if(int.alpha < 0 || int.alpha.err < 0){
    stop("[use_DRAC4flint] Error: data$flint$Dr$alpha or data$flint$Dr$alpha.err is  < 0.")
  }

  if(is.null(int.beta) || is.null(int.beta.err)){
    int.beta <- "X"
    int.beta.err <- "X"
  }else if(!is.numeric(int.beta) || !is.numeric(int.beta.err)){
    stop("[use_DRAC4flint] Error: data$flint$Dr$beta or data$flint$Dr$beta.err is not of type 'numeric'.")
  }else if(int.beta < 0 || int.beta.err < 0){
    stop("[use_DRAC4flint] Error: data$flint$Dr$beta or data$flint$Dr$beta.err is  < 0.")
  }

  if(is.null(int.gamma) || is.null(int.gamma.err)){
    int.gamma <- "X"
    int.gamma.err <- "X"
  }else if(!is.numeric(int.gamma) || !is.numeric(int.gamma.err)){
    stop("[use_DRAC4flint] Error: data$flint$Dr$gamma or data$flint$Dr$gamma.err is not of type 'numeric'.")
  }else if(int.gamma < 0 || int.gamma.err < 0){
    stop("[use_DRAC4flint] Error: data$flint$Dr$gamma or data$flint$Dr$gamma.err is  < 0.")
  }

  # Concentration
  if (is.null(int.U) || is.null(int.Th) || is.null(int.K)){
    int.U <- "X"
    int.U.err <- "X"
    int.Th <- "X"
    int.Th.err <- "X"
    int.K <- "X"
    int.K.err <- "X"
    int.Rb <- "X"
    int.Rb.err <- "X"

    int.K2Rb <- "N"

  }else{
    if(!is.numeric(int.U) || !is.numeric(int.U.err)){
      stop("[use_DRAC4flint] Error: data$flint$Dr$U or data$flint$Dr$U.err is not of type 'numeric'.")
    }else if(int.U < 0 || int.U.err < 0){
      stop("[use_DRAC4flint] Error: data$flint$Dr$U or data$flint$Dr$U.err is  < 0.")
    }

    if(!is.numeric(int.Th) || !is.numeric(int.Th.err)){
      stop("[use_DRAC4flint] Error: data$flint$Dr$Th or data$flint$Dr$Th.err is not of type 'numeric'.")
    }else if(int.Th < 0 || int.Th.err < 0){
      stop("[use_DRAC4flint] Error: data$flint$Dr$Th or data$flint$Dr$Th.err is  < 0.")
    }

    if(!is.numeric(int.K) || !is.numeric(int.K.err)){
      stop("[use_DRAC4flint] Error: data$flint$Dr$K or data$flint$Dr$K.err is not of type 'numeric'.")
    }else if(int.K < 0 || int.K.err < 0){
      stop("[use_DRAC4flint] Error: data$flint$Dr$K or data$flint$Dr$K.err is  < 0.")
    }

    if(is.null(int.Rb) || is.null(int.Rb.err)){
      int.K2Rb <- "Y"
      int.Rb <- "X"
      int.Rb.err <- "X"

    }else if(!is.numeric(int.Rb) || !is.numeric(int.Rb.err)){
      stop("[use_DRAC4flint] Error: data$flint$Dr$Rb or data$flint$Dr$Rb.err is not of type 'numeric'.")
    }else if(int.Rb < 0 || int.Rb.err < 0){
      stop("[use_DRAC4flint] Error: data$flint$Dr$Rb or data$flint$Dr$Rb.err is  < 0.")
    }else if(is.null(int.K2Rb)){
      int.K2Rb <- "N"
    }else if(!is.logical(int.K2Rb)){
      stop("[use_DRAC4flint] Error: 'data$flint$Dr$K2Rb' is not of type 'logical'.")
    }else if(int.K2Rb){
      int.K2Rb <- "Y"
    }else{
      int.K2Rb <- "N"
    }
  }

  ## ------------------------------------------------------------------------------------------

  # External dose rate
  # By concentration
  ext.U <- int.U
  ext.U.err <- int.U.err
  ext.Th <- int.Th
  ext.Th.err <- int.Th.err
  ext.K <- int.K
  ext.K.err <- int.K.err
  ext.Rb <- int.Rb
  ext.Rb.err <- int.Rb.err

  ext.K2Rb <- int.K2Rb


  # Direct evaluation
  ext.alpha <- int.alpha
  ext.alpha.err <- int.alpha.err
  ext.beta <- int.beta
  ext.beta.err <- int.beta.err
  ext.gamma <-int.gamma
  ext.gamma.err <- int.gamma.err

  ## ------------------------------------------------------------------------------------------

  # Environmental dose rate
  # By concentration
  env.U <- data$sediment$Dr$U
  env.U.err <- data$sediment$Dr$U.err
  env.Th <- data$sediment$Dr$Th
  env.Th.err <- data$sediment$Dr$Th.err
  env.K <- data$sediment$Dr$K
  env.K.err <- data$sediment$Dr$K.err
  env.Rb <- data$sediment$Dr$Rb
  env.Rb.err <- data$sediment$Dr$Rb.err

  env.K2Rb <- data$sediment$Dr$K2Rb


  # Direct evaluation
  env.alpha <- data$sediment$Dr$alpha
  env.alpha.err <- data$sediment$Dr$alpha.err
  env.beta <- data$sediment$Dr$beta
  env.beta.err <- data$sediment$Dr$beta.err
  env.gamma <- data$sediment$Dr$gamma
  env.gamma.err <- data$sediment$Dr$gamma.err


  ## ------------------------------------------------------------------------------------------
  # Check values

  # Dose rate direct estimation
  if(is.null(env.alpha) || is.null(env.alpha.err)){
    env.alpha <- "X"
    env.alpha.err <- "X"
  }else if(!is.numeric(env.alpha) || !is.numeric(env.alpha.err)){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$alpha or data$sediment$Dr$alpha.err is not of type 'numeric'.")
  }else if(env.alpha < 0 || env.alpha.err < 0){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$alpha or data$sediment$Dr$alpha.err is  < 0.")
  }

  if(is.null(env.beta) || is.null(env.beta.err)){
    env.beta <- "X"
    env.beta.err <- "X"
  }else if(!is.numeric(env.beta) || !is.numeric(env.beta.err)){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$beta or data$sediment$Dr$beta.err is not of type 'numeric'.")
  }else if(env.beta < 0 || env.beta.err < 0){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$beta or data$sediment$Dr$beta.err is  < 0.")
  }

  if(is.null(env.gamma) || is.null(env.gamma.err)){
    env.gamma <- "X"
    env.gamma.err <- "X"
  }else if(!is.numeric(env.gamma) || !is.numeric(env.gamma.err)){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$gamma or data$sediment$Dr$gamma.err is not of type 'numeric'.")
  }else if(env.gamma < 0 || env.gamma.err < 0){
    stop("[use_DRAC4flint] Error: data$sediment$Dr$gamma or data$sediment$Dr$gamma.err is  < 0.")
  }

  # Dose rate base on concentration
  if (is.null(env.U) || is.null(env.Th) || is.null(env.K)){
    env.U <- "X"
    env.U.err <- "X"
    env.Th <- "X"
    env.Th.err <- "X"
    env.K <- "X"
    env.K.err <- "X"
    env.Rb <- "X"
    env.Rb.err <- "X"

    env.K2Rb <- "N"

  }else{
    if(!is.numeric(env.U) || !is.numeric(env.U.err)){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$U or data$sediment$Dr$U.err is not of type 'numeric'.")
    }else if(env.U < 0 || env.U.err < 0){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$U or data$sediment$Dr$U.err is  < 0.")
    }

    if(!is.numeric(env.Th) || !is.numeric(env.Th.err)){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$Th or data$sediment$Dr$Th.err is not of type 'numeric'.")
    }else if(env.Th < 0 || env.Th.err < 0){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$Th or data$sediment$Dr$Th.err is  < 0.")
    }

    if(!is.numeric(env.K) || !is.numeric(env.K.err)){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$K or data$sediment$Dr$K.err is not of type 'numeric'.")
    }else if(env.K < 0 || env.K.err < 0){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$K or data$sediment$Dr$K.err is  < 0.")
    }

    if(is.null(env.Rb)){
      env.K2Rb <- "Y"
      env.Rb <- "X"
      env.Rb.err <- "X"

    }else if(!is.numeric(env.Rb) || !is.numeric(env.Rb.err)){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$Rb or data$sediment$Dr$Rb.err is not of type 'numeric'.")
    }else if(env.Rb < 0 || env.Rb.err < 0){
      stop("[use_DRAC4flint] Error: data$sediment$Dr$Rb or data$sediment$Dr$Rb.err is  < 0.")
    }else if(is.null(env.K2Rb)){
      env.K2Rb <- "N"
    }else if(!is.logical(env.K2Rb)){
      stop("[use_DRAC4flint] Error: 'data$sediment$Dr$K2Rb' is not of type 'logical'.")
    }else if(env.K2Rb){
      env.K2Rb <- "Y"
    }else{
      env.K2Rb <- "N"
    }
  }
  ## ------------------------------------------------------------------------------------------

  # Flint information

  # Internal information
  int.size.min <- data$flint$info$grain.size.min
  int.size.max <- data$flint$info$grain.size.max

  int.etch.min <- data$flint$info$grain.etch.min
  int.etch.max <- data$flint$info$grain.etch.max

  int.a.value <- data$flint$info$a.value
  int.a.value.err <- data$flint$info$a.value.err

  ## ------------------------------------------------------------------------------------------
  # Check values

  if(is.null(int.size.min) || is.null(int.size.max)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.size.min' or 'data$flint$info$grain.size.max' is null.")
  }else if(!is.numeric(int.size.min) || !is.numeric(int.size.max)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.size.min' or 'data$flint$info$grain.size.max' is not of type 'numeric'.")
  }else if(int.size.min > int.size.max){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.size.min' > 'data$flint$info$grain.size.max'.")
  }else if(int.size.min < 1){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.size.min' < 1 [um].")
  }else if(int.size.max > 1000){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.size.min' > 1000 [um].")
  }

  if(is.null(int.etch.min) || is.null(int.etch.max)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.etch.min' or 'data$flint$info$grain.etch.max' is null.")
  }else if(!is.numeric(int.etch.min) || !is.numeric(int.etch.max)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.etch.min' or 'data$flint$info$grain.etch.max' is not of type 'numeric'.")
  }else if(int.etch.min > int.etch.max){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.etch.min' > 'data$flint$info$grain.etch.max'.")
  }else if(int.etch.min < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.etch.min' < 0 [um].")
  }else if(int.etch.max > 30){
    stop("[use_DRAC4flint] Error: 'data$flint$info$grain.etch.min' > 30 [um].")
  }

  if(is.null(int.a.value) || is.null(int.a.value.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$a.value' or 'data$flint$info$a.value.err' is null.")
  }else if(!is.numeric(int.a.value) || !is.numeric(int.a.value.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$a.value' or 'data$flint$info$a.value.err' is not of type 'numeric'.")
  }else if(int.a.value < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$a.value' < 0 [um].")
  }else if(int.a.value.err < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$a.value.err' < 0 [um].")
  }
  ## ------------------------------------------------------------------------------------------

  # 'External' information
  max.etch <- 30

  ext.etch.min <- max.etch
  ext.etch.max <- max.etch

  ext.water <- data$flint$info$water.content
  ext.water.err <- data$flint$info$water.content.err

  ext.density <- data$flint$info$density
  ext.density.err <- data$flint$info$density.err

  ## ------------------------------------------------------------------------------------------
  # Check values

  if(is.null(ext.water) || is.null(ext.water.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$water.content' or 'data$flint$info$water.content.err' is null.")
  }else if(!is.numeric(ext.water) || !is.numeric(ext.water.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$water.content' or 'data$flint$info$water.content.err' is not of type 'numeric'.")
  }else if(ext.water < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$water.content' < 0 [um].")
  }else if(ext.water.err < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$water.content.err' < 0 [um].")
  }

  if(is.null(ext.density) || is.null(ext.density.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$density' or 'data$flint$info$density.err' is null.")
  }else if(!is.numeric(ext.density) || !is.numeric(ext.density.err)){
    stop("[use_DRAC4flint] Error: 'data$flint$info$density' or 'data$flint$info$density.err' is not of type 'numeric'.")
  }else if(ext.density < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$density' < 0 [um].")
  }else if(ext.density.err < 0){
    stop("[use_DRAC4flint] Error: 'data$flint$info$density.err' < 0 [um].")
  }
  ## ------------------------------------------------------------------------------------------

  # Sediment information

  # Env information
  env.water <- data$sediment$info$water.content
  env.water.err <- data$sediment$info$water.content.err

  env.density <- data$sediment$info$density
  env.density.err <- data$sediment$info$density.err

  env.scale4shallow.depth <- data$sediment$info$scale4shallow.depth

  ## ------------------------------------------------------------------------------------------
  # Check values

  if(is.null(env.water) || is.null(env.water.err)){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$water.content' or 'data$sediment$info$water.content.err' is null.")
  }else if(!is.numeric(env.water) || !is.numeric(env.water.err)){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$water.content' or 'data$sediment$info$water.content.err' is not of type 'numeric'.")
  }else if(env.water < 0){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$water.content' < 0 [um].")
  }else if(env.water.err < 0){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$water.content.err' < 0 [um].")
  }

  if(is.null(env.density) || is.null(env.density.err)){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$density' or 'data$sediment$info$density.err' is null.")
  }else if(!is.numeric(env.density) || !is.numeric(env.density.err)){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$density' or 'data$sediment$info$density.err' is not of type 'numeric'.")
  }else if(env.density < 0){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$density' < 0 [um].")
  }else if(env.density.err < 0){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$density.err' < 0 [um].")
  }

  if(is.null(env.scale4shallow.depth)){
    env.scale4shallow.depth <- 'N'
    warning("[use_DRAC4flint] Error: 'data$sediment$info$scale4shallow.depth' is null.")

  }else if(!is.logical(env.scale4shallow.depth)){
    stop("[use_DRAC4flint] Error: 'data$sediment$info$scale4shallow.depth' is not of type 'logical'.")
  }else if(env.scale4shallow.depth){
    env.scale4shallow.depth <- "Y"
  }else{
    env.scale4shallow.depth <- "N"
  }
  ## ------------------------------------------------------------------------------------------

  # Cosmic dose rate

  #Theoretical

  depth <- data$cosmic$depth
  depth.err <- data$cosmic$depth.err

  latitude <- data$cosmic$latitude
  longitude <- data$cosmic$longitude
  altitude <- data$cosmic$altitude

  #Direct measurement
  cosmic <- data$cosmic$Dr
  cosmic.err <- data$cosmic$Dr.err

  # For calc_CosmicDoseRate
  corr.fieldChanges <- data$cosmic$corr.fieldChanges

  ## ------------------------------------------------------------------------------------------
  # Check values

  if(is.null(depth) || is.null(depth.err)){
    depth <- "X"
    depth.err <- "X"

  }else if(!is.numeric(depth) || !is.numeric(depth.err)){
    stop("[use_DRAC4flint] Error: 'data$cosmic$depth' or 'data$cosmic$depth.err' is not of type 'numeric'.")
  }else if(depth < 0){
    stop("[use_DRAC4flint] Error: 'data$cosmic$depth' < 0 [m].")
  }else if(depth.err < 0){
    stop("[use_DRAC4flint] Error: 'data$cosmic$depth.err' < 0 [m].")
  }

  if(is.null(latitude) || is.null(longitude) || is.null(altitude)){
    latitude <- "X"
    longitude <- "X"
    altitude <- "X"

  }else if(!is.numeric(latitude) || !is.numeric(longitude) || !is.numeric(altitude) ){
    stop("[use_DRAC4flint] Error: 'data$cosmic$latitude', 'data$cosmic$longitude' or 'data$cosmic$altitude' is not of type 'numeric'.")
  }else if(latitude < -90 || latitude > 90){
    stop("[use_DRAC4flint] Error: 'data$cosmic$latitude' has to be between -90\u00b0 and +90\u00b0.")
  }else if(longitude < -180 || longitude > 180){
    stop("[use_DRAC4flint] Error: 'data$cosmic$longitude' has to be between -180\u00b0 and +180\u00b0.")
  }else if(altitude > 5000){
    stop("[use_DRAC4flint] Error: 'data$cosmic$altitude' has to be < 5000 m.")
  }

  if(is.null(cosmic) || is.null(cosmic.err)){
    cosmic <- "X"
    cosmic.err <- "X"
  }else if(!is.numeric(cosmic) || !is.numeric(cosmic.err)){
    stop("[use_DRAC4flint] Error: 'data$cosmic$Dr' or 'data$cosmic$Dr.err' is not of type 'numeric'.")
  }else if(cosmic < 0){
    stop("[use_DRAC4flint] Error: 'data$cosmic$Dr' < 0 [Gy].")
  }else if(cosmic.err < 0){
    stop("[use_DRAC4flint] Error: 'data$cosmic$Dr.err' < 0 [Gy].")
  }

  if(is.null(corr.fieldChanges)){
    corr.fieldChanges <- FALSE
  }else if(!is.logical(corr.fieldChanges)){
    stop("[use_DRAC4flint] Error: 'data$cosmic$corr.fieldChanges' is not of type 'logical'.")
  }

  ## ------------------------------------------------------------------------------------------
  # First run, Internal & External dose (flint)

  int.input <- template_DRAC(notification = notification)

  int.input$`Project ID` <- project
  int.input$`Sample ID` <- sample
  int.input$Mineral <- mineral
  int.input$`Conversion factors` <- conversion.factors

  int.input$`Internal U (ppm)` <-  int.U
  int.input$`errInternal U (ppm)` <- int.U.err
  int.input$`Internal Th (ppm)` <- int.Th
  int.input$`errInternal Th (ppm)` <- int.Th.err
  int.input$`Internal K (%)` <- int.K
  int.input$`errInternal K (%)` <- int.K.err
  int.input$`Rb (ppm)` <- int.Rb
  int.input$`errRb (ppm)` <- int.Rb.err
  int.input$`Calculate internal Rb from K conc?`<- int.K2Rb

  int.input$`ExternalU (ppm)` <-  ext.U
  int.input$`errExternal U (ppm)` <- ext.U.err
  int.input$`External Th (ppm)` <- ext.Th
  int.input$`errExternal Th (ppm)` <- ext.Th.err
  int.input$`External K (%)` <- ext.K
  int.input$`errExternal K (%)` <- ext.K.err
  int.input$`External Rb (ppm)` <- ext.Rb
  int.input$`errExternal Rb (ppm)` <- ext.Rb.err
  int.input$`Calculate external Rb from K conc?` <- ext.K2Rb

  int.input$`Grain size min (microns)` <- int.size.min
  int.input$`Grain size max (microns)` <- int.size.max

  int.input$`alpha-Grain size attenuation` <- int.alpha.attenuation
  int.input$`beta-Grain size attenuation ` <- int.beta.attenuation

  int.input$`Etch depth min (microns)` <- int.etch.min
  int.input$`Etch depth max (microns)` <- int.etch.max

  int.input$`beta-Etch depth attenuation factor` <- int.etch.beta.attenuation

  int.input$`a-value` <- int.a.value
  int.input$`erra-value` <- int.a.value.err

  int.input$`Water content ((wet weight - dry weight)/dry weight) %` <- ext.water
  int.input$`errWater content %` <- ext.water.err

  int.input$`Depth (m)` <- depth
  int.input$`errDepth (m)` <- depth.err

  int.input$`Overburden density (g cm-3)` <- env.density
  int.input$`errOverburden density (g cm-3)` <- env.density.err

  int.input$`Latitude (decimal degrees)` <- latitude
  int.input$`Longitude (decimal degrees)` <- longitude
  int.input$`Altitude (m)` <- altitude

  int.input$`De (Gy)` <- De
  int.input$`errDe (Gy)` <- De.err

  int.output <- use_DRAC(int.input,verbose=notification)

  # --------------------------------------------------------------------------------------------------------------------------

  # Second run, Evironmental dose (sediment)

  env.input <- template_DRAC(notification = FALSE)

  env.input$`Project ID` <- project
  env.input$`Sample ID` <- sample
  env.input$Mineral <- mineral
  env.input$`Conversion factors` <- conversion.factors

  env.input$`ExternalU (ppm)` <-  env.U
  env.input$`errExternal U (ppm)` <- env.U.err
  env.input$`External Th (ppm)` <- env.Th
  env.input$`errExternal Th (ppm)` <- env.Th.err
  env.input$`External K (%)` <- env.K
  env.input$`errExternal K (%)` <- env.K.err
  int.input$`External Rb (ppm)` <- env.Rb
  int.input$`errExternal Rb (ppm)` <- env.Rb.err
  env.input$`Calculate external Rb from K conc?` <- env.K2Rb

  env.input$`Scale gammadoserate at shallow depths?` <- env.scale4shallow.depth

  env.input$`Grain size min (microns)` <- int.size.min
  env.input$`Grain size max (microns)` <- int.size.max

  env.input$`a-value` <- int.a.value
  env.input$`erra-value` <- int.a.value.err

  env.input$`Water content ((wet weight - dry weight)/dry weight) %` <- env.water
  env.input$`errWater content %` <- env.water.err

  env.input$`Depth (m)` <- depth
  env.input$`errDepth (m)` <- depth.err
  env.input$`Overburden density (g cm-3)` <- env.density
  env.input$`errOverburden density (g cm-3)` <- env.density.err
  env.input$`Latitude (decimal degrees)` <- latitude
  env.input$`Longitude (decimal degrees)` <- longitude
  env.input$`Altitude (m)` <- altitude

  env.input$`De (Gy)` <- De
  env.input$`errDe (Gy)` <- De.err

  env.output <- use_DRAC(env.input,verbose=FALSE)

  # --------------------------------------------------------------------------------------------------------------------------

  # Combining the 2 runs...

  temp.int.alpha <- as.numeric(int.output$DRAC$highlights$`Internal Dry alphadoserate (Gy.ka-1)`[1])
  temp.int.alpha.err <- as.numeric(int.output$DRAC$highlights$`Internal Dry erralphadoserate (Gy.ka-1)`[1])

  temp.int.beta <- as.numeric(int.output$DRAC$highlights$`Internal Dry betadoserate (Gy.ka-1)`[1])
  temp.int.beta.err <- as.numeric(int.output$DRAC$highlights$`Internal Dry errbetadoserate (Gy.ka-1)`[1])

  temp.ext.alpha <- as.numeric(int.output$DRAC$highlights$`Water corrected alphadoserate`[1])
  temp.ext.alpha.err <- as.numeric(int.output$DRAC$highlights$`Water corrected erralphadoserate`[1])

  temp.ext.beta <- as.numeric(int.output$DRAC$highlights$`Water corrected betadoserate`[1])
  temp.ext.beta.err <- as.numeric(int.output$DRAC$highlights$`Water corrected errbetadoserate`[1])

  temp.ext.gamma <- as.numeric(int.output$DRAC$highlights$`Water corrected gammadoserate (Gy.ka-1)`[1])
  temp.ext.gamma.err <- as.numeric(int.output$DRAC$highlights$`Water corrected errgammadoserate (Gy.ka-1)`[1])

  temp.env.alpha <- as.numeric(env.output$DRAC$highlights$`Water corrected alphadoserate`[1])
  temp.env.alpha.err <- as.numeric(env.output$DRAC$highlights$`Water corrected erralphadoserate`[1])

  temp.env.beta <-  as.numeric(env.output$DRAC$highlights$`Water corrected betadoserate`[1])
  temp.env.beta.err <- as.numeric(env.output$DRAC$highlights$`Water corrected errbetadoserate`[1])

  temp.env.gamma <- as.numeric(env.output$DRAC$highlights$`Water corrected gammadoserate (Gy.ka-1)`[1])
  temp.env.gamma.err <- as.numeric(env.output$DRAC$highlights$`Water corrected errgammadoserate (Gy.ka-1)`[1])

  temp.cosmic <- as.numeric(env.output$DRAC$highlights$`Cosmicdoserate (Gy.ka-1)`[1])
  temp.cosmic.err <- as.numeric(env.output$DRAC$highlights$`errCosmicdoserate (Gy.ka-1)`[1])

  # --------------------------------------------------------------------------------------------------------------------------

  # Estimation of external gamma over/under.estimation

  if(temp.env.gamma > temp.ext.gamma){
    gamma.over <- abs(temp.env.gamma-temp.ext.gamma)/temp.ext.gamma*100

    message.gamma <- paste("The gamma dose rate coming from the flint is probably overestimated (+",round(gamma.over,2)," %).",sep="")
  }else if(temp.env.gamma > temp.ext.gamma){
    gamma.under <- abs(temp.env.gamma-temp.ext.gamma)/temp.ext.gamma

    message.gamma <- paste("The gamma dose rate coming from the flint is probably underestimated (-",round(gamma.under,2)," %).",sep="")
  }else{
    message.gamma <- "The environment and the flint produce similar gamma dose rate. \n"
  }

  comment <- c(comment, gamma=message.gamma)

  # DRAC results
  # Partial dose rate
  DRAC.int.Dr <- temp.int.alpha+temp.int.beta
  DRAC.int.Dr.err <- sqrt(sum(temp.int.alpha.err^2, temp.int.beta.err^2))

  DRAC.ext.Dr <- temp.ext.alpha+temp.ext.beta
  DRAC.ext.Dr.err <- sqrt(sum(temp.ext.alpha.err^2, temp.ext.beta.err^2))

  DRAC.env.Dr <- temp.env.gamma+temp.cosmic
  DRAC.env.Dr.err <- sqrt(sum(temp.env.gamma.err^2, temp.cosmic.err^2))

  DRAC.alpha.Dr <- temp.int.alpha+temp.ext.alpha
  DRAC.alpha.Dr.err <- sqrt(sum(temp.int.alpha.err^2, temp.ext.alpha.err^2))

  DRAC.beta.Dr <- temp.int.beta+temp.ext.beta
  DRAC.beta.Dr.err <- sqrt(sum(temp.int.beta.err^2, temp.ext.beta.err^2))

  DRAC.gamma.Dr <- temp.env.gamma
  DRAC.gamma.Dr.err <- temp.env.gamma.err

  DRAC.cosmic.Dr <- temp.cosmic
  DRAC.cosmic.Dr.err <- temp.cosmic.err

  # Total dose rate
  DRAC.Dr <- sum(DRAC.alpha.Dr,
                 DRAC.beta.Dr,
                 DRAC.gamma.Dr,
                 DRAC.cosmic.Dr)

  DRAC.Dr.err <- sqrt(sum(DRAC.alpha.Dr.err^2,
                          DRAC.beta.Dr.err^2,
                          DRAC.gamma.Dr.err^2,
                          DRAC.cosmic.Dr.err^2))

  # Age
  DRAC.age <- De/DRAC.Dr
  DRAC.age.err <- DRAC.age*sqrt(sum((De.err/De)^2,(DRAC.Dr.err/DRAC.Dr)^2))


  # output
  DRAC.result <- list(Age = DRAC.age,
                      Age.err =DRAC.age.err,
                      De=De,
                      De.err = De.err,
                      Dr = DRAC.Dr,
                      Dr.err = DRAC.Dr.err,
                      int.Dr = DRAC.int.Dr,
                      int.Dr.err = DRAC.int.Dr.err,
                      ext.Dr = DRAC.ext.Dr,
                      ext.Dr.err = DRAC.ext.Dr.err,
                      env.Dr = DRAC.env.Dr,
                      env.Dr.err = DRAC.env.Dr.err,
                      alpha.Dr = DRAC.alpha.Dr,
                      alpha.Dr.err = DRAC.alpha.Dr.err,
                      beta.Dr = DRAC.beta.Dr,
                      beta.Dr.err = DRAC.beta.Dr.err,
                      gamma.Dr = DRAC.gamma.Dr,
                      gamma.Dr.err = DRAC.gamma.Dr.err,
                      cosmic.Dr = DRAC.cosmic.Dr,
                      cosmic.Dr.err = DRAC.cosmic.Dr.err)

  # Correction of the cosmic dose rate using 'calc_CosmicDoseRate'

  R.cosmic <- calc_CosmicDoseRate(depth = depth,
                                  density = env.density,
                                  latitude = latitude,
                                  longitude = longitude,
                                  altitude = altitude,
                                  est.age = DRAC.age,
                                  error = DRAC.cosmic.Dr.err/DRAC.cosmic.Dr,
                                  corr.fieldChanges = corr.fieldChanges)

  R.cosmic.Dr <- R.cosmic$summary$dc
  R.cosmic.Dr.err <- R.cosmic$summary$dc*R.cosmic$args$error

  # R results

  # Partial dose rate
  R.int.Dr <- temp.int.alpha+temp.int.beta
  R.int.Dr.err <- sqrt(sum(temp.int.alpha.err^2, temp.int.beta^2))

  R.ext.Dr <- temp.ext.alpha+temp.ext.beta
  R.ext.Dr.err <- sqrt(sum(temp.ext.alpha.err^2, temp.ext.beta^2))

  R.env.Dr <- temp.env.gamma+R.cosmic.Dr
  R.env.Dr.err <- sqrt(sum(temp.env.gamma.err^2, R.cosmic.Dr.err^2))

  R.alpha.Dr <- temp.int.alpha+temp.ext.alpha
  R.alpha.Dr.err <- sqrt(sum(temp.int.alpha.err^2, temp.ext.alpha^2))

  R.beta.Dr <- temp.int.beta+temp.ext.beta
  R.beta.Dr.err <- sqrt(sum(temp.int.beta.err^2, temp.ext.beta^2))

  R.gamma.Dr <- temp.env.gamma
  R.gamma.Dr.err <- temp.env.gamma.err

  # Total dose rate
  R.Dr <- sum(R.alpha.Dr,
              R.beta.Dr,
              R.gamma.Dr,
              R.cosmic.Dr)

  R.Dr.err <- sqrt(sum(R.alpha.Dr.err^2,
                       R.beta.Dr.err^2,
                       R.gamma.Dr.err^2,
                       R.cosmic.Dr.err^2))

  # Age
  R.age <- De/R.Dr
  R.age.err <- R.age*sqrt(sum((De.err/De)^2,(R.Dr.err/R.Dr)^2))



  # Output
  R.result <- list(Age = R.age,
                   Age.err =R.age.err,
                   De=De,
                   De.err = De.err,
                   Dr = R.Dr,
                   Dr.err = R.Dr.err,
                   int.Dr = R.int.Dr,
                   int.Dr.err = R.int.Dr.err,
                   ext.Dr = R.ext.Dr,
                   ext.Dr.err = R.ext.Dr.err,
                   env.Dr = R.env.Dr,
                   env.Dr.err = R.env.Dr.err,
                   alpha.Dr = R.alpha.Dr,
                   alpha.Dr.err = R.alpha.Dr.err,
                   beta.Dr = R.beta.Dr,
                   beta.Dr.err = R.beta.Dr.err,
                   gamma.Dr = R.gamma.Dr,
                   gamma.Dr.err = R.gamma.Dr.err,
                   cosmic.Dr = R.cosmic.Dr,
                   cosmic.Dr.err = R.cosmic.Dr.err)


  age.CE <- date - R.age*1000
  age.CE.err <- R.age.err*1000

  if(age.CE > 0){
    message.CE <- paste("The heating of the flint occured around", round(age.CE), "\u00b1", round(age.CE.err), "CE.")
  }else if(abs(age.CE)<100){
    message.CE <- paste("The heating of the flint occured around", round(abs(age.CE),-1), "\u00b1", round(age.CE.err,-1), " BCE.")
  }else if(abs(age.CE)<1000){
    message.CE <- paste("The heating of the flint occured around", round(abs(age.CE),-2), "\u00b1", round(age.CE.err,-2), " BCE.")
  }else{
    message.CE <- paste("The heating of the flint occured around", round(abs(age.CE),-3), "\u00b1", round(age.CE.err,-3), " BCE.")
  }

  message.project <- paste("For the sample", sample, "of the project", project)
  message.De <- paste("The equivalent dose is:", round(De,3), "\u00b1", round(De.err,3), "Gy.")
  message.Dr <- paste("The dose rate is:", round(R.Dr,3), "\u00b1", round(R.Dr.err,3), "Gy/ka.")
  message.Age <- paste("The age is estimated as:", round(R.age,3), "\u00b1", round(R.age.err,3), "ka.")

  output <- paste("\t [use_DRAC4flint] \n ",
                  "\t \n",
                  paste("\t", message.project, "\n"),
                  "\t --------------------------------------------------------- \n ",
                  paste("\t", message.De, "\n"),
                  paste("\t", message.Dr, "\n"),
                  paste("\t", message.gamma, "\n"),
                  "\t --------------------------------------------------------- \n",
                  paste("\t", message.Age,"\n"),
                  paste("\n \t", message.CE, "\n"))


  cat(output)

  result <- list(age=R.age,
                 age.err=R.age.err,
                 Dr=R.Dr,
                 Dr.err=R.Dr.err,
                 DRAC = DRAC.result,
                 R = R.result,
                 comment=comment)

  new.TLum.Results.use_DRAC4flint <- set_TLum.Results(data = result)

  return (new.TLum.Results.use_DRAC4flint)

  return(result)
}