calc_Huntley2006.Rd

% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/calc_Huntley2006.R
\name{calc_Huntley2006}
\alias{calc_Huntley2006}
\title{Apply the Huntley (2006) model}
\usage{
calc_Huntley2006(
  data,
  LnTn = NULL,
  rhop,
  ddot,
  readerDdot,
  normalise = TRUE,
  fit.method = c("EXP", "GOK"),
  lower.bounds = c(-Inf, -Inf, -Inf, -Inf),
  summary = TRUE,
  plot = TRUE,
  ...
)
}
\arguments{
\item{data}{\link{data.frame} (\strong{required}):
A \code{data.frame} with one of the following structures:
\itemize{
\item A \strong{three column} data frame with numeric values on a) dose (s), b) \code{LxTx} and
c) \code{LxTx} error.
\item If a \strong{two column} data frame is provided it is automatically
assumed that errors on \code{LxTx} are missing. A third column will be attached
with an arbitrary 5 \\% error on the provided \code{LxTx} values.
\item Can also be a \strong{wide table}, i.e. a \link{data.frame} with a number of columns divisible by 3
and where each triplet has the aforementioned column structure.
}

\if{html}{\out{<div class="sourceCode">}}\preformatted{                        (optional)
     | dose (s)| LxTx | LxTx error |
     |  [ ,1]  | [ ,2]|    [ ,3]   |
     |---------|------|------------|
[1, ]|  0      | LnTn | LnTn error | (optional, see arg 'LnTn')
[2, ]|  R1     | L1T1 | L1T1 error |
 ... |    ...  |  ... |     ...    |
[x, ]|  Rx     | LxTx | LxTx error |

}\if{html}{\out{</div>}}

\strong{NOTE:} The function assumes the first row of the function to be the
\code{Ln/Tn}-value. If you want to provide more than one \code{Ln/Tn}-value consider
using the argument \code{LnTn}.}

\item{LnTn}{\link{data.frame} (\strong{optional}):
This argument should \strong{only} be used to provide more than one \code{Ln/Tn}-value.
It assumes a two column data frame with the following structure:

\if{html}{\out{<div class="sourceCode">}}\preformatted{     |  LnTn  |  LnTn error  |
     |  [ ,1] |      [ ,2]   |
     |--------|--------------|
[1, ]| LnTn_1 | LnTn_1 error |
[2, ]| LnTn_2 | LnTn_2 error |
 ... |   ...  |      ...     |
[x, ]| LnTn_x | LnTn_x error |
}\if{html}{\out{</div>}}

The function will calculate a \strong{mean} \code{Ln/Tn}-value and uses either the
standard deviation or the highest individual error, whichever is larger. If
another mean value (e.g. a weighted mean or median) or error is preferred,
this value must be calculated beforehand and used in the first row in the
data frame for argument \code{data}.

\strong{NOTE:} If you provide \code{LnTn}-values with this argument the data frame
for the \code{data}-argument \strong{must not} contain any \code{LnTn}-values!}

\item{rhop}{\link{numeric} (\strong{required}):
The density of recombination centres (\eqn{\rho}') and its error (see Huntley 2006),
given as numeric vector of length two. Note that \eqn{\rho}' must \strong{not} be
provided as the common logarithm. Example: \code{rhop = c(2.92e-06, 4.93e-07)}.}

\item{ddot}{\link{numeric} (\strong{required}):
Environmental dose rate and its error, given as a numeric vector of length two.
Expected unit: Gy/ka. Example: \code{ddot = c(3.7, 0.4)}.}

\item{readerDdot}{\link{numeric} (\strong{required}):
Dose rate of the irradiation source of the OSL reader and its error,
given as a numeric vector of length two.
Expected unit: Gy/s. Example: \code{readerDdot = c(0.08, 0.01)}.}

\item{normalise}{\link{logical} (\emph{with default}): If \code{TRUE} (the default) all measured and computed \eqn{\frac{L_x}{T_x}} values are normalised by the pre-exponential factor \code{A} (see details).}

\item{fit.method}{\link{character} (\emph{with default}):
Fit function of the dose response curve. Can either be \code{EXP} (the default)
or \code{GOK}. Note that \code{EXP} (single saturating exponential) is the original
function the model after Huntley (2006) and Kars et al. (2008) was
designed to use. The use of a general-order kinetics function (\code{GOK})
is an experimental adaptation of the model and should be used
with great care.}

\item{lower.bounds}{\link{numeric} (\emph{with default}):
A vector of length 4 that contains the lower bound values to be applied
when fitting the models with \link[minpack.lm:nlsLM]{minpack.lm::nlsLM}. In most cases, the
default values (\code{c(-Inf, -Inf, -Inf, -Inf)}) are appropriate for finding
a best fit, but sometimes it may be useful to restrict the lower bounds to
e.g. \code{c(0, 0, 0, 0)}. The values of the vectors are, respectively, for
parameters \code{a}, \code{D0}, \code{c} and \code{d} in that order (parameter \code{d} is ignored
when \code{fit.method = "EXP"}). More details can be found in
\link{fit_DoseResponseCurve}.}

\item{summary}{\link{logical} (\emph{with default}):
If \code{TRUE} (the default) various parameters provided by the user
and calculated by the model are added as text on the right-hand side of the
plot.}

\item{plot}{\link{logical} (\emph{with default}):
enables/disables plot output.}

\item{...}{Further parameters:
\itemize{
\item \code{verbose} \link{logical}: Show or hide console output
\item \code{n.MC} \link{numeric}: Number of Monte Carlo iterations (default = \code{10000}).
\strong{Note} that it is generally advised to have a large number of Monte Carlo
iterations for the results to converge. Decreasing the number of iterations
will often result in unstable estimates.
}

All other arguments are passed to \link{plot} and \link{fit_DoseResponseCurve} (in particular
\code{mode} for the fit mode and \code{fit.force_through_origin})}
}
\value{
An \linkS4class{RLum.Results} object is returned:

Slot: \strong{@data}\cr

\tabular{lll}{
\strong{OBJECT} \tab \strong{TYPE} \tab \strong{COMMENT}\cr
\code{results} \tab \link{data.frame} \tab results of the of Kars et al. 2008 model \cr
\code{data} \tab \link{data.frame} \tab original input data \cr
\code{Ln} \tab \link{numeric} \tab Ln and its error \cr
\code{LxTx_tables} \tab \code{list} \tab A \code{list} of \code{data.frames} containing data on dose,
LxTx and LxTx error for each of the dose response curves.
Note that these \strong{do not} contain the natural Ln signal, which is provided separately. \cr
\code{fits} \tab \code{list} \tab A \code{list} of \code{nls} objects produced by \link[minpack.lm:nlsLM]{minpack.lm::nlsLM} when fitting the dose response curves \cr
}

Slot: \strong{@info}\cr

\tabular{lll}{
\strong{OBJECT} \tab \strong{TYPE} \tab \strong{COMMENT} \cr
\code{call} \tab \code{call} \tab the original function call \cr
\code{args} \tab \code{list} \tab arguments of the original function call \cr
}
}
\description{
A function to calculate the expected sample specific fraction of saturation
based on the model of Huntley (2006) using the approach as implemented
in Kars et al. (2008) or Guralnik et al. (2015).
}
\details{
This function applies the approach described in Kars et al. (2008) or Guralnik et al. (2015),
which are both developed from the model of Huntley (2006) to calculate the expected sample
specific fraction of saturation of a feldspar and also to calculate fading
corrected age using this model. \eqn{\rho}' (\code{rhop}), the density of recombination
centres, is a crucial parameter of this model and must be determined
separately from a fading measurement. The function \link{analyse_FadingMeasurement}
can be used to calculate the sample specific \eqn{\rho}' value.

\strong{Kars et al. (2008) - Single saturating exponential}

To apply the approach after Kars et al. (2008) use \code{fit.method = "EXP"}.

Firstly, the unfaded \eqn{D_0} value is determined through applying equation 5 of
Kars et al. (2008) to the measured \eqn{\frac{L_x}{T_x}} data as a function of irradiation
time, and fitting the data with a single saturating exponential of the form:

\deqn{\frac{L_x}{T_x}(t^*) = A  \phi(t^*) \{1 - exp(-\frac{t^*}{D_0}))\}}

where

\deqn{\phi(t^*) = exp(-\rho' ln(1.8  \tilde{s} t^*)^3)}

after King et al. (2016) where \eqn{A} is a pre-exponential factor,
\eqn{t^*} (s) is the irradiation time, starting at the mid-point of
irradiation (Auclair et al. 2003) and \eqn{\tilde{s}} (\eqn{3\times10^{15}} s\eqn{^{-1}}) is the athermal frequency factor after Huntley (2006). \cr

Using fit parameters \eqn{A} and \eqn{D_0}, the function then computes a natural dose
response curve using the environmental dose rate, \eqn{\dot{D}} (Gy/s) and equations
\verb{[1]} and \verb{[2]}. Computed \eqn{\frac{L_x}{T_x}} values are then fitted using the
\link{fit_DoseResponseCurve} function and the laboratory measured LnTn can then
be interpolated onto this curve to determine the fading corrected
\eqn{D_e} value, from which the fading corrected age is calculated.

\strong{Guralnik et al. (2015) - General-order kinetics}

To apply the approach after Guralnik et al. (2015) use \code{fit.method = "GOK"}.

The approach of Guralnik et al. (2015) is very similar to that of
Kars et al. (2008), but instead of using a single saturating exponential
the model fits a general-order kinetics function of the form:

\deqn{\frac{L_x}{T_x}(t^*) = A \phi (t^*)(1 - (1 + (\frac{1}{D_0}) t^* c)^{-1/c})}

where \eqn{A}, \eqn{\phi}, \eqn{t^*} and \eqn{D_0} are the same as above and \eqn{c} is a
dimensionless kinetic order modifier (cf. equation 10 in
Guralnik et al., 2015).

\strong{Level of saturation}

The \link{calc_Huntley2006} function also calculates the level of saturation (\eqn{\frac{n}{N}})
and the field saturation (i.e. athermal steady state, (n/N)_SS) value for
the sample under investigation using the sample specific \eqn{\rho}',
unfaded \eqn{D_0} and \eqn{\dot{D}} values, following the approach of Kars et al. (2008).

The computation is done using 1000 equally-spaced points in the interval
[0.01, 3]. This can be controlled by setting option \code{rprime}, such as
in \code{rprime = seq(0.01, 3, length.out = 1000)} (the default).

\strong{Uncertainties}

Uncertainties are reported at \eqn{1\sigma} and are assumed to be normally
distributed and are estimated using Monte-Carlo re-sampling (\code{n.MC = 1000})
of \eqn{\rho}' and \eqn{\frac{L_x}{T_x}} during dose response curve fitting, and of \eqn{\rho}'
in the derivation of (\eqn{n/N}) and (n/N)_SS.

\strong{Age calculated from 2D0 of the simulated natural DRC}

In addition to the age calculated from the equivalent dose derived from
\eqn{\frac{L_n}{T_n}} projected on the simulated natural dose response curve (DRC), this function
also calculates an age from twice the characteristic saturation dose (\code{D0})
of the simulated natural DRC. This can be a useful information for
(over)saturated samples (i.e., no intersect of \eqn{\frac{L_n}{T_n}} on the natural DRC)
to obtain at least a "minimum age" estimate of the sample. In the console
output this value is denoted by \emph{"Age @2D0 (ka):"}.
}
\note{
This function has BETA status, in particular for the GOK implementation. Please verify
your results carefully
}
\section{Function version}{
 0.4.5
}

\examples{

## Load example data (sample UNIL/NB123, see ?ExampleData.Fading)
data("ExampleData.Fading", envir = environment())

## (1) Set all relevant parameters
# a. fading measurement data (IR50)
fading_data <- ExampleData.Fading$fading.data$IR50

# b. Dose response curve data
data <- ExampleData.Fading$equivalentDose.data$IR50

## (2) Define required function parameters
ddot <- c(7.00, 0.004)
readerDdot <- c(0.134, 0.0067)

# Analyse fading measurement and get an estimate of rho'.
# Note that the RLum.Results object can be directly used for further processing.
# The number of MC runs is reduced for this example
rhop <- analyse_FadingMeasurement(fading_data, plot = TRUE, verbose = FALSE, n.MC = 10)

## (3) Apply the Kars et al. (2008) model to the data
kars <- calc_Huntley2006(
 data = data,
 rhop = rhop,
 ddot = ddot,
 readerDdot = readerDdot,
 n.MC = 25)


\dontrun{
# You can also provide LnTn values separately via the 'LnTn' argument.
# Note, however, that the data frame for 'data' must then NOT contain
# a LnTn value. See argument descriptions!
LnTn <- data.frame(
 LnTn = c(1.84833, 2.24833),
 nTn.error = c(0.17, 0.22))

LxTx <- data[2:nrow(data), ]

kars <- calc_Huntley2006(
 data = LxTx,
 LnTn = LnTn,
 rhop = rhop,
 ddot = ddot,
 readerDdot = readerDdot,
 n.MC = 25)
}
} 

\section{How to cite}{
King, G.E., Burow, C., Kreutzer, S., 2024. calc_Huntley2006(): Apply the Huntley (2006) model. Function version 0.4.5. In: Kreutzer, S., Burow, C., Dietze, M., Fuchs, M.C., Schmidt, C., Fischer, M., Friedrich, J., Mercier, N., Philippe, A., Riedesel, S., Autzen, M., Mittelstrass, D., Gray, H.J., Galharret, J., Colombo, M., 2024. Luminescence: Comprehensive Luminescence Dating Data Analysis. R package version 0.9.25. https://r-lum.github.io/Luminescence/
}

\references{
Kars, R.H., Wallinga, J., Cohen, K.M., 2008. A new approach towards anomalous fading correction for feldspar
IRSL dating-tests on samples in field saturation. Radiation Measurements 43, 786-790. doi:10.1016/j.radmeas.2008.01.021

Guralnik, B., Li, B., Jain, M., Chen, R., Paris, R.B., Murray, A.S., Li, S.-H., Pagonis, P.,
Herman, F., 2015. Radiation-induced growth and isothermal decay of infrared-stimulated luminescence
from feldspar. Radiation Measurements 81, 224-231.

Huntley, D.J., 2006. An explanation of the power-law decay of luminescence.
Journal of Physics: Condensed Matter 18, 1359-1365. doi:10.1088/0953-8984/18/4/020

King, G.E., Herman, F., Lambert, R., Valla, P.G., Guralnik, B., 2016.
Multi-OSL-thermochronometry of feldspar. Quaternary Geochronology 33, 76-87. doi:10.1016/j.quageo.2016.01.004

\strong{Further reading}

Morthekai, P., Jain, M., Cunha, P.P., Azevedo, J.M., Singhvi, A.K., 2011. An attempt to correct
for the fading in million year old basaltic rocks. Geochronometria 38(3), 223-230.
}
\author{
Georgina E. King, University of Lausanne (Switzerland) \cr
Christoph Burow, University of Cologne (Germany) \cr
Sebastian Kreutzer, Ruprecht-Karl University of Heidelberg (Germany)
, RLum Developer Team}
\keyword{datagen}