TLdating: Tools for Thermoluminescence Dating

#' MAAD protocol for TL dating
#'
#' Function to estimate the ED in TL dating using the MAAD protocol. \cr
#' It provides an estimation of the palaeodose (Q) and/or the supralinearity correction (I).
#' The equivalent dose (ED) is estimated by the addition of Q and I. \cr
#' See details for more information.
#'
#'
#' @param object
#'  \code{\linkS4class{TLum.Analysis}} (\bold{required}): object containing the TL curves used for the ED estimation.
#' @param eval.Tmin
#'  \link{integer} (\bold{required}): Temperature ( °C) of the lower boundary for the signal integration.
#' @param eval.Tmax
#'  \link{integer} (\bold{required}): Temperature ( °C) of the upper boundary for the signal integration.
#' @param rejection.criteria
#'  \link{list} (with default): list containing the rejection criteria (in \%). See details.
#' @param fitting.parameters
#'  \link{list} (with default): list containing the fitting parameters. See details.
#' @param plotting.parameters
#'  \link{list} (with default): list containing the plotting parameters. See details.
#'
#' @details
#' This function estimates the equivalent dose for the thermoluminescence dating with the MAAD protocol.
#' It can provide an estimation of the palaeodose (Q) and the supralinearity correction (I) simultaniously or separately.
#' These are estimated using the growth curve approach (QC) (Aitken, 1985) and the dose plateau approach (DP).
#' Both approaches should provide a similar result. The equivalent dose is estimated by the addition of Q and I\cr
#' The Lx/Tx matrix is estimated using \link{calc_TL.LxTx}. \cr
#' The average TL curves for each dose step are estimate using \link{calc_TL.MAAD.average}. \cr
#' The plateau test values are estimated using \link{calc_TL.plateau}. \cr
#'
#' \bold{Rejection criteria} \cr
#' The rejection criteria are: \cr
#' \describe{
#'  \item{\code{testdose.error}}{
#'    \link{numeric}: Maximum error accepted on Tx (in \%).}
#'  \item{\code{paleodose.error}}{
#'    \link{numeric}: Maximum error accepted on Lx (in \%).}
#'  }
#'
#' \bold{Fitting parameters} \cr
#' The fitting parameters are:  \cr
#' \describe{
#'  \item{\code{method}}{
#'    \link{character}: Fitting method (\code{LIN}, \code{EXP}, \code{EXP+LIN} or \code{EXP+EXP}).}
#'  \item{\code{fit.weighted}}{
#'    \link{logical}: If the fitting is weighted or not.}
#'  \item{\code{fit.use.slope}}{
#'    \link{logical}: If the slope of the Q growth curve is reused for the supralinearity correction.}
#'  \item{\code{fit.aDoses.min}}{
#'    \link{numeric}: Lowest additive dose used for the fitting.}
#'  \item{\code{fit.aDoses.max}}{
#'    \link{numeric}: Highest additive dose used for the fitting.}
#'  \item{\code{fit.rDoses.min}}{
#'    \link{numeric}: Lowest regenerative dose used for the fitting.}
#'  \item{\code{fit.rDoses.max}}{
#'    \link{numeric}: Highest regenerative dose used for the fitting.}
#' }
#' See also \link{calc_TL.MAAD.fit.Q} and \link{calc_TL.MAAD.fit.I}. \cr
#'
#' \bold{Plotting parameters} \cr
#' The plotting parameters are:  \cr
#' \describe{
#'  \item{\code{plot.Tmin}}{
#'    \link{numeric}: Lowest temperature plotted.}
#'  \item{\code{plot.Tmax}}{
#'    \link{numeric}: Highest temperature plotted.}
#'  \item{\code{no.plot}}{
#'    \link{logical}: If \code{TRUE}, the results will not be plotted.}
#' }
#' See also \link{plot_TL.MAAD}. \cr
#'
#' @return
#'  The results are plotted using \link{plot_TL.MAAD}. \cr
#'
#'  The function also provides a \linkS4class{TLum.Results} containing: \cr
#'  \describe{
#'    \item{\code{De.GC}}{
#'      \link{list}: Results obtained with the dose plateau approach and their uncertainties
#'      (\code{De}, \code{De.error}, \code{Q}, \code{Q.error}, \code{I}, \code{I.error})}
#'    \item{\code{De.DP}}{
#'      \link{list}: Results obtained with the growth curve approach and their uncertainties
#'      (\code{De}, \code{De.error}, \code{Q}, \code{Q.error}, \code{I}, \code{I.error})}
#'    \item{\code{LnLxTnTx.table}}{
#'      \link{matrix}: Lx/Tx values}
#'    \item{\code{RC.Status}}{
#'    \link{character}: The acceptance result. }
#'  }
#'
#'
#' @references
#'  Aitken, M.J. (1985) Thermoluminescence Dating, Academic Press, London \cr
#'
#' @seealso
#'  \link{calc_TL.LxTx},
#'  \link{calc_TL.plateau},
#'  \link{calc_TL.MAAD.average},
#'  \link{calc_TL.MAAD.separate},
#'  \link{calc_TL.MAAD.fit.I},
#'  \link{calc_TL.MAAD.fit.Q},
#'  \link{analyse_TL.SAR}.
#'
#'
#' @author David Strebler, University of Cologne (Germany)
#'
#' @examples
#'
#' ##load data
#'
#' ##perform analysis
#'
#' @export analyse_TL.MAAD


analyse_TL.MAAD <- function(
  object,

  eval.Tmin,

  eval.Tmax,

  rejection.criteria = list(testdose.error = 10,
                            paleodose.error = 10),

  fitting.parameters=list(fit.method="LIN",
                          fit.weighted=FALSE,
                          fit.use.slope=FALSE,
                          fit.aDoses.min=0,
                          fit.aDoses.max=NA,
                          fit.rDoses.min=0,
                          fit.rDoses.max=NA),

  plotting.parameters=list(plot.Tmin=0,
                           plot.Tmax=NA,
                           no.plot=FALSE)

  ){

  # ------------------------------------------------------------------------------
  # Integrity Check
  # ------------------------------------------------------------------------------

  if (missing(object)){
    stop("[analyse_TL.SAR] Error: Input 'object' is missing.")

  }else if (!is(object,"TLum.Analysis")){
    stop("[analyse_TL.MAAD] Error: Input 'object' is not of type 'TLum.Analysis'.")
  }

  if (missing(eval.Tmin)){
    stop("[analyse_TL.MAAD] Error: Input 'eval.Tmin' is missing.")

  }else if(!is.numeric(eval.Tmin)){
    stop("[calc_TL.MAAD.average] Error: Input 'eval.Tmin' is not of type 'numeric'.")
  }

  if (missing(eval.Tmax)){
    stop("[analyse_TL.MAAD] Error: Input 'eval.Tmax' is missing.")

  }else if(!is.numeric(eval.Tmax)){
    stop("[calc_TL.MAAD.average] Error: Input 'eval.Tmax' is not of type 'numeric'.")
  }
  # ------------------------------------------------------------------------------

  sample.name <- as.character(object@records[[1]]@metadata$SAMPLE)

  data <- calc_TL.LxTx(object)

  temperatures  <- get_TLum.Results(data, "Temperatures")

  nPoints <- length(temperatures)
  Tmax <- max(temperatures)
  Tstep <- Tmax/nPoints

  # ------------------------------------------------------------------------------
  # Value check

  if(eval.Tmin < 0){
    stop("[analyse_TL.MAAD] Error: eval.Tmin < 0.")
  }

  if(eval.Tmax > Tmax){
    stop(paste("[analyse_TL.MAAD] Error: eval.Tmax >", Tmax, "(Tmax)."))
  }

  if(eval.Tmin > eval.Tmax){
    stop("[analyse_TL.MAAD] Error: eval.Tmin > eval.Tmax.")
  }

  if(eval.Tmin < 0){
    stop("[analyse_TL.MAAD] Error: eval.Tmin < 0.")
  }
  # ------------------------------------------------------------------------------

  eval.min <- ceiling(eval.Tmin/Tstep)
  eval.max <- floor(eval.Tmax/Tstep)

  dTypes <- get_TLum.Results(data, "Datatypes")

  doses <- get_TLum.Results(data, "Doses")

  Lx <- as.matrix(get_TLum.Results(data, "Lx"))
  Lx.error <- as.matrix(get_TLum.Results(data, "Lx.error"))

  Tx <- as.matrix(get_TLum.Results(data, "Tx"))
  Tx.error <- as.matrix(get_TLum.Results(data, "Tx.error"))

  LxTx <- as.matrix(get_TLum.Results(data, "LxTx"))
  LxTx.error <- as.matrix(get_TLum.Results(data, "LxTx.error"))


  # ------------------------------------------------------------------------------
  # Separation of aLx and rLx
  # ------------------------------------------------------------------------------

  temp.data <- calc_TL.MAAD.separate(doses=doses,
                                     Lx=Lx,
                                     Lx.error=Lx.error,
                                     dTypes=dTypes)

  aDoses.f <- get_TLum.Results(temp.data, "aDoses")
  aNames.f <- get_TLum.Results(temp.data, "aNames")

  rDoses.f <- get_TLum.Results(temp.data, "rDoses")
  rNames.f <- get_TLum.Results(temp.data, "rNames")

  aLx.f <- as.matrix(get_TLum.Results(temp.data, "aLx"))
  aLx.f.error <- as.matrix(get_TLum.Results(temp.data, "aLx.error"))

  rLx.f <- as.matrix(get_TLum.Results(temp.data, "rLx"))
  rLx.f.error <- as.matrix(get_TLum.Results(temp.data, "rLx.error"))

  temp.data <- calc_TL.MAAD.separate(doses=doses,
                                     Lx=Tx,
                                     Lx.error=Tx.error,
                                     dTypes=dTypes)

  aTx.f <- as.matrix(get_TLum.Results(temp.data, "aLx"))
  aTx.f.error <- as.matrix(get_TLum.Results(temp.data, "aLx.error"))

  rTx.f <- as.matrix(get_TLum.Results(temp.data, "rLx"))
  rTx.f.error <- as.matrix(get_TLum.Results(temp.data, "rLx.error"))

  temp.data <- calc_TL.MAAD.separate(doses=doses,
                                     Lx=LxTx,
                                     Lx.error=LxTx.error,
                                     dTypes=dTypes)

  aLxTx.f <- as.matrix(get_TLum.Results(temp.data, "aLx"))
  aLxTx.f.error <- as.matrix(get_TLum.Results(temp.data, "aLx.error"))

  rLxTx.f <- as.matrix(get_TLum.Results(temp.data, "rLx"))
  rLxTx.f.error <- as.matrix(get_TLum.Results(temp.data, "rLx.error"))


  # ------------------------------------------------------------------------------
  # Average for aLx and rLx
  # ------------------------------------------------------------------------------

  # Additive doses
  aDoses.a <- vector()
  aNames.a <- vector()

  aLx.a <- vector()
  aLx.a.error <- vector()

  aTx.a <- vector()
  aTx.a.error <- vector()

  aLxTx.a  <- vector()
  aLxTx.a.error <- vector()

  if(length(aLx.f) > 0){
    temp.data <- calc_TL.MAAD.average(names=aNames.f,
                                      doses=aDoses.f,
                                      Lx=aLx.f,
                                      Lx.error=aLx.f.error)

    aDoses.a <- get_TLum.Results(temp.data, "doses")
    aNames.a <- get_TLum.Results(temp.data, "names")

    aLx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
    aLx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))

    if(length(aTx.f) > 0){


      temp.data <- calc_TL.MAAD.average(names=aNames.f,
                                        doses=aDoses.f,
                                        Lx=aTx.f,
                                        Lx.error=aTx.f.error)

      aTx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
      aTx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))
    }

    temp.data <- calc_TL.MAAD.average(names=aNames.f,
                                      doses=aDoses.f,
                                      Lx=aLxTx.f,
                                      Lx.error=aLxTx.f.error)

    aLxTx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
    aLxTx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))
  }

  #Regenerative doses
  rDoses.a <- vector()
  rNames.a <- vector()

  rLx.a <- vector()
  rLx.a.error <- vector()

  rTx.a <- vector()
  rTx.a.error <- vector()

  rLxTx.a  <- vector()
  rLxTx.a.error <- vector()

  if(length(rLx.f) > 0){
    temp.data <- calc_TL.MAAD.average(names=rNames.f,
                                    doses=rDoses.f,
                                    Lx=rLx.f,
                                    Lx.error=rLx.f.error)

    rDoses.a <- get_TLum.Results(temp.data, "doses")
    rNames.a <- get_TLum.Results(temp.data, "names")

    rLx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
    rLx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))

    if(length(rTx.f) > 0){
      temp.data <- calc_TL.MAAD.average(names=rNames.f,
                                        doses=rDoses.f,
                                        Lx=rTx.f,
                                        Lx.error=rTx.f.error)

      rTx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
      rTx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))
    }

    temp.data <- calc_TL.MAAD.average(names=rNames.f,
                                      doses=rDoses.f,
                                      Lx=rLxTx.f,
                                      Lx.error=rLxTx.f.error)

    rLxTx.a <- as.matrix(get_TLum.Results(temp.data, "Lx"))
    rLxTx.a.error <- as.matrix(get_TLum.Results(temp.data, "Lx.error"))
  }

  # ------------------------------------------------------------------------------
  # Plateau test
  # ------------------------------------------------------------------------------

  #Plateau test (Additive step)
    #aLx
  if(length(aLx.a) > 0){
    aLx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aLx.a[,aNames.a == "N"],
                                                      Ln.error=aLx.a.error[,aNames.a == "N"],
                                                      Lx=aLx.a[,aNames.a != "N" & aDoses.a != 0],
                                                      Lx.error=aLx.a.error[,aNames.a != "N"& aDoses.a != 0]),
                                      "LnLx")
  }
    #atx
  if(length(aTx.a) > 0){
    aTx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aTx.a[,aNames.a == "N"],
                                                      Ln.error=aTx.a.error[,aNames.a == "N"],
                                                      Lx=aTx.a[,aNames.a != "N" & aDoses.a != 0],
                                                      Lx.error=aTx.a.error[,aNames.a != "N"& aDoses.a != 0]),
                                      "LnLx")
  }
    #aLxTx
  if(length(aLxTx.a) > 0){
    aLxTx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aLxTx.a[,aNames.a == "N"],
                                                        Ln.error=aLxTx.a.error[,aNames.a == "N"],
                                                        Lx=aLxTx.a[,aNames.a != "N" & aDoses.a != 0],
                                                        Lx.error=aLxTx.a.error[,aNames.a != "N"& aDoses.a != 0]),
                                        "LnLx")
  }

  #Plateau test (Additive step)
    #rLx
  if(length(rLx.a) > 0){
    if(length(aLx.a)>0){
      rLx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aLx.a[,aNames.a == "N"],
                                                        Ln.error=aLx.a.error[,aNames.a == "N"],
                                                        Lx=rLx.a[,rNames.a != "N" & rDoses.a != 0],
                                                        Lx.error=rLx.a.error[,rNames.a != "N"& rDoses.a != 0]),
                                        "LnLx")
    }
    else{
      rLx.a.plateau <- rLx.a
    }
  }
    #rTx
  if(length(rTx.a) > 0){
    if(length(aTx.a)>0){
      rTx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aTx.a[,aNames.a == "N"],
                                                        Ln.error=aTx.a.error[,aNames.a == "N"],
                                                        Lx=rTx.a[,rNames.a != "N" & rDoses.a != 0],
                                                        Lx.error=rTx.a.error[,rNames.a != "N"& rDoses.a != 0]),
                                        "LnLx")
    }else{
      rTx.a.plateau<-rTx.a
    }
  }
    #rLxTx
  if(length(rLxTx.a) > 0){
    if(length(aLxTx.a)>0){
      rLxTx.a.plateau <- get_TLum.Results(calc_TL.plateau(Ln=aLxTx.a[,aNames.a == "N"],
                                                          Ln.error=aLxTx.a.error[,aNames.a == "N"],
                                                          Lx=rLxTx.a[,rNames.a != "N" & rDoses.a != 0],
                                                          Lx.error=rLxTx.a.error[,rNames.a != "N"& rDoses.a != 0]),
                                          "LnLx")
    }else{
      rLxTx.a.plateau<-rLxTx.a
    }
  }

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

  # ------------------------------------------------------------------------------
  # Selection criteria
  # ------------------------------------------------------------------------------

  fit.aDoses.min <- fitting.parameters$fit.aDoses.min
  fit.aDoses.max <- fitting.parameters$fit.aDoses.max

  fit.rDoses.min <- fitting.parameters$fit.rDoses.min
  fit.rDoses.max <- fitting.parameters$fit.rDoses.max

  paleodose.error <- rejection.criteria$paleodose.error
  testdose.error <- rejection.criteria$testdose.error

  # ------------------------------------------------------------------------------
  # Value check
  # additive doses
  if(length(aDoses.a) > 0){
    if(is.null(fit.aDoses.min) || is.na(fit.aDoses.min) || fit.aDoses.min < min(aDoses.a)){
      fit.aDoses.min <- min(aDoses.a[aDoses.a != 0])
    }

    if(is.null(fit.aDoses.max) || is.na(fit.aDoses.max) || fit.aDoses.max > max(aDoses.a)){
      fit.aDoses.max <- max(aDoses.a)
    }

    if(fit.aDoses.max < fit.aDoses.min){
      stop("[analyse_TL.SAR] Error: fit.aDoses.max < fit.aDoses.min")
    }
  }else{
    fit.aDoses.min <- NULL
    fit.aDoses.max <- NULL
  }

  # regenerative doses
  if(length(rDoses.a) > 0){
    if(is.null(fit.rDoses.min) || is.na(fit.rDoses.min) || fit.rDoses.min < min(rDoses.a)){
      fit.rDoses.min <- min(rDoses.a[rDoses.a != 0])
    }

    if(is.null(fit.rDoses.max) || is.na(fit.rDoses.max) || fit.rDoses.max > max(rDoses.a)){
      fit.rDoses.max <- max(rDoses.a)
    }

    if(fit.rDoses.max < fit.rDoses.min){
      stop("[analyse_TL.SAR] Error: fit.rDoses.max < fit.rDoses.min")
    }
  }else{
    fit.rDoses.min <- NULL
    fit.rDoses.max <- NULL
  }

  # rejection.criteria
  if(is.null(paleodose.error) || !is.finite(paleodose.error)){
    stop("[analyse_TL.SAR] Error: paleodose.error is missing.")
  }else if(paleodose.error<0){
    paleodose.error <- abs(paleodose.error)
    warning("[analyse_TL.SAR] Warning: paleodose.error is negative.")
  }

  if(is.null(testdose.error) || !is.finite(testdose.error)){
    stop("[analyse_TL.SAR] Error: paleodose.error is missing.")
  }else if(testdose.error<0){
    testdose.error <- abs(testdose.error)
    warning("[analyse_TL.SAR] Warning: testdose.error is negative.")
  }
  #---------------------------------------------------------
  Lx.error.criteria <- round(abs(paleodose.error/100),digits=3)
  Tx.error.criteria <- round(abs(testdose.error/100),digits=3)


  #Max paleodose error (addivite curve)
  if(length(aLx.f) > 0){
    temp.aLx.lim <- aLx.f[eval.min:eval.max, aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max]
    temp.aLx.lim.error <- aLx.f.error[eval.min:eval.max, aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max]
    temp.aLx.lim.error.r <- abs(temp.aLx.lim.error/temp.aLx.lim)

    temp.aLx.lim.error.r[!is.finite(temp.aLx.lim.error.r)] <- NA

    aLx.error.r.GC <- vector()

    for(i in 1:ncol(temp.aLx.lim.error.r)){
      temp.aLx.error.r.GC <- mean(temp.aLx.lim.error.r[,i],na.rm = TRUE)

      aLx.error.r.GC <- c(aLx.error.r.GC, temp.aLx.error.r.GC)
    }

    aLx.error.max <- max(aLx.error.r.GC, na.rm = TRUE)
  }else{
    aLx.error.max <- NA
  }

  #Max testdose error (additive curve)
  if(length(aTx.f) > 0){
    temp.aTx.lim <- aTx.f[eval.min:eval.max, aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max]
    temp.aTx.lim.error <- aTx.f.error[eval.min:eval.max, aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max]
    temp.aTx.lim.error.r <- abs(temp.aTx.lim.error/temp.aTx.lim)

    temp.aTx.lim.error.r[!is.finite(temp.aTx.lim.error.r)] <- NA

    aTx.error.r.GC <- vector()

    for(i in 1:ncol(temp.aTx.lim.error.r)){
      temp.aTx.error.r.GC <- mean(temp.aTx.lim.error.r[,i],na.rm = TRUE)

      aTx.error.r.GC <- c(aTx.error.r.GC, temp.aTx.error.r.GC)
    }

    aTx.error.max <- max(aTx.error.r.GC, na.rm = TRUE)
  }else{
    aTx.error.max <- NA
  }

  #Max paleodose error (regenerative curve)
  if(length(rLx.f) > 0){
    temp.rLx.lim <- rLx.f[eval.min:eval.max, rDoses.f >= fit.rDoses.min & rDoses.f <= fit.rDoses.max & rDoses.f != 0]
    temp.rLx.lim.error <- rLx.f.error[eval.min:eval.max, rDoses.f >= fit.rDoses.min & rDoses.f <= fit.rDoses.max & rDoses.f != 0]
    temp.rLx.lim.error.r <- abs(temp.rLx.lim.error/temp.rLx.lim)

    temp.rLx.lim.error.r[!is.finite(temp.rLx.lim.error.r)] <- NA

    rLx.error.r.GC <- vector()

    for(i in 1:ncol(temp.rLx.lim.error.r)){
      temp.rLx.error.r.GC <- mean(temp.rLx.lim.error.r[,i],na.rm = TRUE)

      rLx.error.r.GC <- c(rLx.error.r.GC, temp.rLx.error.r.GC)
    }

    rLx.error.max <- max(rLx.error.r.GC, na.rm = TRUE)
  }else{
    rLx.error.max <- NA
  }

  #Max testdose error (regenerative curve)
  if(length(rTx.f) > 0){
    temp.rTx.lim <- rTx.f[eval.min:eval.max,rDoses.a >= fit.rDoses.min & rDoses.a <= fit.rDoses.max]
    temp.rTx.lim.error <- rTx.f.error[eval.min:eval.max,rDoses.a >= fit.rDoses.min & rDoses.a <= fit.rDoses.max ]
    temp.rTx.lim.error.r <- abs(temp.rTx.lim.error/temp.rTx.lim)

    temp.rTx.lim.error.r[!is.finite(temp.rTx.lim.error.r)] <- NA

    rTx.error.r.GC <- vector()

    for(i in 1:ncol(temp.rTx.lim.error.r)){
      temp.rTx.error.r.GC <- mean(temp.rTx.lim.error.r[,i],na.rm = TRUE)

      rTx.error.r.GC <- c(rTx.error.r.GC, temp.rTx.error.r.GC)
    }

    rTx.error.max <- max(rTx.error.r.GC, na.rm = TRUE)
  }else{
    rTx.error.max <- NA
  }

  # Max paleodose error test (additive dose)
  if(length(aLx.f) > 0){
    if(aLx.error.max > Lx.error.criteria){
      test.aLx.error <- FALSE
    }else{
      test.aLx.error <- TRUE
    }
  }else{
    test.aLx.error <- TRUE
  }

  # Max testdose error test (additive dose)
  if(length(aTx.f) > 0){
    if(aTx.error.max > Tx.error.criteria){
      test.aTx.error <- FALSE
    }else{
      test.aTx.error <- TRUE
    }
  }else{
    test.aTx.error <- TRUE
  }

  # Max paleodose test (regenerative curve)
  if(length(rLx.f) > 0){
    if(rLx.error.max > Lx.error.criteria){
      test.rLx.error <- FALSE
    }else{
      test.rLx.error <- TRUE
    }
  }else{
    test.rLx.error <- TRUE
  }

  # Max testdose test (regenerative curve)
  if(length(rTx.f) > 0){
    if(rTx.error.max > Tx.error.criteria){
      test.rTx.error <- FALSE
    }else{
      test.rTx.error <- TRUE
    }
  }else{
    test.rTx.error <- TRUE
  }

  # Acceptance result
  if(test.aLx.error && test.aTx.error && test.rLx.error && test.rTx.error) {
    acceptance.result <- "OK"
  }else{
    acceptance.result <- "FAILED"
  }

  rejection.values <-list(aLx.error.max=aLx.error.max,
                          aTx.error.max=aTx.error.max,
                          rLx.error.max=rLx.error.max,
                          rTx.error.max=rTx.error.max,
                          test.aLx.error=test.aLx.error,
                          test.aTx.error=test.aTx.error,
                          test.rLx.error=test.rLx.error,
                          test.rTx.error=test.rTx.error)

  #----------------------------------------------------------------------------------------------------------------
  # D_e estimation
  #----------------------------------------------------------------------------------------------------------------

  # ------------------------------------------------------------------------------
  # Palaeodose - Q
  # ------------------------------------------------------------------------------

  if(length(aLx.a) > 0){

    # Growth curve (GC)
    # aLxTx.a.lim <- aLxTx.a[eval.min:eval.max,]
    # aLxTx.a.error.lim <- aLxTx.a.error[eval.min:eval.max,]
    aLxTx.f.lim <- aLxTx.f[eval.min:eval.max,]
    aLxTx.f.error.lim <- aLxTx.f.error[eval.min:eval.max,]

    GC.aLxTx <- vector()
    GC.aLxTx.error <- vector()
    GC.aDoses <- aDoses.f
    GC.aNames <- aNames.f

    # for(i in 1:length(aDoses.a)){
    for(i in 1:length(aDoses.f)){

      # temp.LxTx <- aLxTx.a.lim[,i]
      # temp.LxTx.error <- aLxTx.a.error.lim[,i]
      temp.LxTx <- aLxTx.f.lim[,i]
      temp.LxTx.error <- aLxTx.f.error.lim[,i]
      temp.LxTx.w <- 1/(temp.LxTx.error^2)

      GC.aLxTx[i] <- sum(temp.LxTx.w*temp.LxTx,na.rm=TRUE)/sum(temp.LxTx.w,na.rm=TRUE)
      GC.aLxTx.error[i] <- 1/sqrt(sum(temp.LxTx.w,na.rm = TRUE))
    }

    # Q (GC)

      #Data
    # temp.bool <- aDoses.a >= fit.aDoses.min & aDoses.a <= fit.aDoses.max
    temp.bool <- aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max

    # temp.doses <- aDoses.a[temp.bool]
    temp.doses <- aDoses.f[temp.bool]

    temp.LxTx <- GC.aLxTx[temp.bool]
    temp.LxTx.error <- GC.aLxTx.error[temp.bool]

    temp.fit <- calc_TL.MAAD.fit.Q(LxTx = temp.LxTx,
                                LxTx.error = temp.LxTx.error,
                                doses = temp.doses,
                                fitting.parameters=fitting.parameters)

    GC.Q <- get_TLum.Results(temp.fit, "GC")

    Q.GC <- get_TLum.Results(temp.fit, "Q")
    Q.GC.error <- get_TLum.Results(temp.fit, "Q.error")
    Q.GC.slope <- get_TLum.Results(temp.fit, "summary")

    # Q (DP)
    Q.DP <- vector()
    Q.DP.error <- vector()
    Q.DP.slope <- list()

    for(i in 1:length(temperatures)){
      # temp.aLxTx <- aLxTx.a[i,]
      # temp.aLxTx.error <- aLxTx.a.error[i,]
      temp.aLxTx <- aLxTx.f[i,]
      temp.aLxTx.error <- aLxTx.f.error[i,]

      # Data
      # temp.bool <- aDoses.a >= fit.aDoses.min & aDoses.a <= fit.aDoses.max
      temp.bool <- aDoses.f >= fit.aDoses.min & aDoses.f <= fit.aDoses.max

      # temp.doses <- aDoses.a[temp.bool]
      temp.doses <- aDoses.f[temp.bool]

      temp.LxTx <- temp.aLxTx[temp.bool]
      temp.LxTx.error <- temp.aLxTx.error[temp.bool]

      # Regression
      temp.fit <- calc_TL.MAAD.fit.Q(LxTx = temp.LxTx,
                                   LxTx.error = temp.LxTx.error,
                                   doses = temp.doses,
                                   fitting.parameters=fitting.parameters)

      temp.GC <- get_TLum.Results(temp.fit, "GC")

      temp.Q <- get_TLum.Results(temp.fit, "Q")
      temp.Q.error <- get_TLum.Results(temp.fit, "Q.error")
      temp.slope <- get_TLum.Results(temp.fit, "summary")

      #save
      Q.DP <- c(Q.DP, temp.Q)
      Q.DP.error <- c(Q.DP.error, temp.Q.error)
      Q.DP.slope[[i]] <- temp.slope
    }

    # Q.DP average
    Q.DP.w <- 1/(Q.DP.error^2)

    Q.DP.lim <- Q.DP[eval.min:eval.max]
    Q.DP.lim.error <- Q.DP.error[eval.min:eval.max]
    Q.DP.lim.w <- 1/(Q.DP.lim.error^2)

    Q.DP.a <- sum(Q.DP.lim.w*Q.DP.lim,na.rm=TRUE)/sum(Q.DP.lim.w,na.rm=TRUE)
    Q.DP.a.error <- 1/sqrt(sum(Q.DP.lim.w, na.rm = TRUE))
  }else{
    GC.aLxTx <- vector()
    GC.aLxTx.error <- vector()

    GC.Q <- vector()

    Q.DP <- vector()
    Q.DP.error <- vector()
    Q.DP.a <- vector()
    Q.DP.a.error <- vector()
    Q.GC <- vector()
    Q.GC.error <- vector()

    Q.DP.slope <- list()
    Q.GC.slope <- list()

    }

  if(length(Q.DP) == 0){
    Q.DP.a <- 0
    Q.DP.a.error <- 0
    Q.GC <- 0
    Q.GC.error <- 0
  }

  # ------------------------------------------------------------------------------
  # supralinearity correction - I
  # ------------------------------------------------------------------------------

  if(length(rLxTx.a) > 0){

    # Growth curve
    # rLxTx.a.w <- 1/(rLxTx.a.error^2)
    #
    # rLxTx.a.lim <- rLxTx.a[eval.min:eval.max,]
    # rLxTx.a.w.lim <- rLxTx.a.w[eval.min:eval.max,]

    rLxTx.f.lim <- rLxTx.f[eval.min:eval.max,]
    rLxTx.f.error.lim <- rLxTx.f.error[eval.min:eval.max,]

    GC.rLxTx <- vector()
    GC.rLxTx.error <- vector()
    GC.rDoses <- rDoses.f
    GC.rNames <- rNames.f

    # for(i in 1:length(rDoses.a)){
    for(i in 1:length(rDoses.f)){

      # temp.LxTx <- rLxTx.a.lim[,i]
      # temp.LxTx.w <- rLxTx.a.w.lim[,i]
      temp.LxTx <- rLxTx.f.lim[,i]
      temp.LxTx.error <- rLxTx.f.error.lim[,i]
      temp.LxTx.w <- 1/(temp.LxTx.error^2)

      GC.rLxTx[i] <- sum(temp.LxTx.w*temp.LxTx, na.rm=TRUE)/sum(temp.LxTx.w, na.rm=TRUE)
      GC.rLxTx.error[i] <- 1/sqrt(sum(temp.LxTx.w,na.rm = TRUE))
    }

    #I (GC)
      # Data
    # temp.bool <- rDoses.a >= fit.rDoses.min & rDoses.a <= fit.rDoses.max
    temp.bool <- rDoses.f >= fit.rDoses.min & rDoses.f <= fit.rDoses.max

    # temp.doses <- rDoses.a[temp.bool]
    temp.doses <- rDoses.f[temp.bool]

    temp.LxTx <- GC.rLxTx[temp.bool]
    temp.LxTx.error <- GC.rLxTx.error[temp.bool]


    if(length(GC.Q)>0){
      temp.slope <- Q.GC.slope

      temp.fit <- calc_TL.MAAD.fit.I(LxTx = temp.LxTx,
                                   LxTx.error = temp.LxTx.error,
                                   doses = temp.doses,
                                   slope = temp.slope,
                                   fitting.parameters=fitting.parameters)
    }else{

      fitting.parameters$fit.use.slope =FALSE
      temp.fit <- calc_TL.MAAD.fit.I(LxTx = temp.LxTx,
                                   LxTx.error = temp.LxTx.error,
                                   doses = temp.doses,
                                   fitting.parameters=fitting.parameters)
    }


    GC.I <- get_TLum.Results(temp.fit, "GC")

    I.GC <- get_TLum.Results(temp.fit, "I")
    I.GC.error <- get_TLum.Results(temp.fit, "I.error")
    I.GC.slope <- get_TLum.Results(temp.fit, "summary")


    # I (DP)
    I.DP <- vector()
    I.DP.error <- vector()
    I.DP.slope <- list()

    for(i in 1:length(temperatures)){
      # temp.rLxTx <- rLxTx.a[i,]
      # temp.rLxTx.error <- rLxTx.a.error[i,]
      temp.rLxTx <- rLxTx.f[i,]
      temp.rLxTx.error <- rLxTx.f.error[i,]


        # selection of the doses used
      # temp.bool <- rDoses.a >= fit.rDoses.min & rDoses.a <= fit.rDoses.max
      temp.bool <- rDoses.f >= fit.rDoses.min & rDoses.f <= fit.rDoses.max

      # temp.dose <- rDoses.a[temp.bool]
      temp.dose <- rDoses.f[temp.bool]

      temp.LxTx <- temp.rLxTx[temp.bool]
      temp.LxTx.error <- temp.rLxTx.error[temp.bool]

      #Regression
      if(length(GC.Q)>0){
        temp.slope <- Q.DP.slope[[i]]

        temp.fit <- calc_TL.MAAD.fit.I(LxTx = temp.LxTx,
                                     LxTx.error = temp.LxTx.error,
                                     doses = temp.dose,
                                     slope = temp.slope,
                                     fitting.parameters=fitting.parameters)
      }else{

        fitting.parameters$fit.use.slope =FALSE
        temp.fit <- calc_TL.MAAD.fit.I(LxTx = temp.LxTx,
                                     LxTx.error = temp.LxTx.error,
                                     doses = temp.dose,
                                     fitting.parameters=fitting.parameters)
      }

      temp.GC <- get_TLum.Results(temp.fit, "GC")

      temp.I <- get_TLum.Results(temp.fit, "I")
      temp.I.error <- get_TLum.Results(temp.fit, "I.error")
      temp.slope <- get_TLum.Results(temp.fit, "summary")

      I.DP <- c(I.DP, temp.I)
      I.DP.error <- c(I.DP.error, temp.I.error)
      I.DP.slope[[i]] <- temp.slope
    }

    # I average
    I.DP.w <- 1/(I.DP.error^2)

    I.DP.lim <- I.DP[eval.min:eval.max]
    I.DP.lim.error <- I.DP.error[eval.min:eval.max]
    I.DP.lim.w <- I.DP.w[eval.min:eval.max]

    I.DP.a <- sum(I.DP.lim.w*I.DP.lim, na.rm=TRUE)/sum(I.DP.lim.w, na.rm=TRUE)
    I.DP.a.error <- 1/sqrt(sum(I.DP.lim.w, na.rm = TRUE))

  }else{
    GC.rLxTx <- vector()
    GC.rLxTx.error <- vector()

    I.DP <- vector()
    I.DP.error <- vector()
    I.DP.a <- vector()
    I.DP.a.error <- vector()
    I.GC <- vector()
    I.GC.error <- vector()
    I.DP.slope <- list()
    I.GC.slope <- list()
  }

  if(length(I.DP) == 0){
    I.DP.a <- 0
    I.DP.a.error <- 0
    I.GC <- 0
    I.GC.error <- 0
  }

  # ------------------------------------------------------------------------------
  # Equivalent dose - De (Q+I)
  # ------------------------------------------------------------------------------
  #De.GC
  De.GC <- sum(Q.GC, I.GC, na.rm = TRUE)
  De.GC.error <- sqrt(sum(Q.GC.error^2,I.GC.error^2, na.rm = TRUE))

  if(!is.finite(De.GC)){
    De.GC <- 0
  }

  if(!is.finite(De.GC.error)){
    De.GC.error <- De.GC
  }

  # De.DP
  De.DP <- sum(Q.DP.a, I.DP.a ,na.rm = TRUE)
  De.DP.error <- sqrt(sum(Q.DP.a.error^2,I.DP.a.error^2, na.rm = TRUE))


  if(!is.finite(De.DP)){
    De.DP <- 0
  }

  if(!is.finite(De.DP.error)){
    De.DP.error <- De.DP
  }

  #----------------------------------------------------------------------------------------------------------------
  #Generate Results
  #----------------------------------------------------------------------------------------------------------------

  new.originator <- as.character(match.call()[[1]])

  new.De.DP <- list(De = De.DP,
                    De.error=De.DP.error,
                    Q = Q.DP.a,
                    Q.error = Q.DP.a.error,
                    I = I.DP.a,
                    I.error = I.DP.a.error)

  new.De.GC <- list(De = De.GC,
                    De.error=De.GC.error,
                    Q = Q.GC,
                    Q.error = Q.GC.error,
                    I = I.GC,
                    I.error = I.GC.error)

  new.data <- list(DP = new.De.DP,
                   GC = new.De.GC,
                   LnLxTnTx.table = LxTx,
                   RC.Status = acceptance.result)

  new.plotData <- list(sample.name=sample.name,
                       fitting.parameters=fitting.parameters,
                       temperatures=temperatures,
                       eval.Tmin=eval.Tmin,
                       eval.Tmax=eval.Tmax,
                       aNames=aNames.a,
                       aDoses=aDoses.a,
                       aLx=aLx.a,
                       aTx=aTx.a,
                       aLxTx=aLxTx.a,
                       aLx.plateau=aLx.a.plateau,
                       aTx.plateau=aTx.a.plateau,
                       aLxTx.plateau=aLxTx.a.plateau,
                       rNames=rNames.a,
                       rDoses=rDoses.a,
                       rLx=rLx.a,
                       rTx=rTx.a,
                       rLxTx=rLxTx.a,
                       rLx.plateau=rLx.a.plateau,
                       rTx.plateau=rTx.a.plateau,
                       rLxTx.plateau=rLxTx.a.plateau,
                       DP.Q.line=Q.DP,
                       DP.Q.line.error=Q.DP.error,
                       GC.Q.slope=Q.GC.slope,
                       GC.Q.line=GC.Q,
                       GC.Q.LxTx=GC.aLxTx,
                       GC.Q.LxTx.error=GC.aLxTx.error,
                       GC.Q.doses = GC.aDoses,
                       GC.Q.names = GC.aNames,
                       DP.I.line=I.DP,
                       DP.I.line.error=I.DP.error,
                       GC.I.slope=I.GC.slope,
                       GC.I.line=GC.I,
                       GC.I.LxTx=GC.rLxTx,
                       GC.I.LxTx.error=GC.rLxTx.error,
                       GC.I.doses=GC.rDoses,
                       GC.I.names=GC.rNames,
                       Q.DP=Q.DP.a,
                       Q.DP.error=Q.DP.a.error,
                       Q.GC=Q.GC,
                       Q.GC.error=Q.GC.error,
                       I.DP=I.DP.a,
                       I.DP.error=I.DP.a.error,
                       I.GC=I.GC,
                       I.GC.error=I.GC.error,
                       De.GC=De.GC,
                       De.GC.error=De.GC.error,
                       De.DP=De.DP,
                       De.DP.error=De.DP.error,
                       rejection.values=rejection.values,
                       plotting.parameters=plotting.parameters
  )

  new.TLum.Results.analyse_TL.MAAD <- set_TLum.Results(originator = new.originator,
                                                       data = new.data,
                                                       plotData =new.plotData)


  #----------------------------------------------------------------------------------------------------------------
  #Plot results
  #----------------------------------------------------------------------------------------------------------------

  no.plot <- plotting.parameters$no.plot

  # ------------------------------------------------------------------------------
  #Check values
  if(is.null(no.plot) || is.na(no.plot) || !is.logical(no.plot)){
    no.plot <- FALSE
  }
  # ------------------------------------------------------------------------------

  if(!no.plot){
    do.call(what = plot_TL.MAAD,
            args = new.plotData)
  }

  #----------------------------------------------------------------------------------------------------------------
  #Export Results
  #----------------------------------------------------------------------------------------------------------------

  return(new.TLum.Results.analyse_TL.MAAD)
}