#' Estimate the amount of grains on an aliquot
#'
#' Estimate the number of grains on an aliquot. Alternatively, the packing
#' density of an aliquot is computed.
#'
#' This function can be used to either estimate the number of grains on an
#' aliquot or to compute the packing density depending on the the arguments
#' provided.
#'
#' The following function is used to estimate the number of grains `n`:
#'
#' \deqn{n = (\pi*x^2)/(\pi*y^2)*d}
#'
#' where `x` is the radius of the aliquot size (microns), `y` is the mean
#' radius of the mineral grains (mm) and `d` is the packing density
#' (value between 0 and 1).
#'
#' **Packing density**
#'
#' The default value for `packing.density` is 0.65, which is the mean of
#' empirical values determined by Heer et al. (2012) and unpublished data from
#' the Cologne luminescence laboratory. If `packing.density = "Inf"` a maximum
#' density of \eqn{\pi/\sqrt12 = 0.9068\ldots} is used. However, note that
#' this value is not appropriate as the standard preparation procedure of
#' aliquots resembles a PECC (*"Packing Equal Circles in a Circle"*) problem
#' where the maximum packing density is asymptotic to about 0.87.
#'
#' **Monte Carlo simulation**
#'
#' The number of grains on an aliquot can be estimated by Monte Carlo simulation
#' when setting `MC = TRUE`. Each of the parameters necessary to calculate
#' `n` (`x`, `y`, `d`) are assumed to be normally distributed with means
#' \eqn{\mu_x, \mu_y, \mu_d} and standard deviations \eqn{\sigma_x, \sigma_y, \sigma_d}.
#'
#' For the mean grain size random samples are taken first from
#' \eqn{N(\mu_y, \sigma_y)}, where \eqn{\mu_y = mean.grain.size} and
#' \eqn{\sigma_y = (max.grain.size-min.grain.size)/4} so that 95\% of all
#' grains are within the provided the grain size range. This effectively takes
#' into account that after sieving the sample there is still a small chance of
#' having grains smaller or larger than the used mesh sizes. For each random
#' sample the mean grain size is calculated, from which random subsamples are
#' drawn for the Monte Carlo simulation.
#'
#' The packing density is assumed
#' to be normally distributed with an empirically determined \eqn{\mu = 0.65}
#' (or provided value) and \eqn{\sigma = 0.18}. The normal distribution is
#' truncated at `d = 0.87` as this is approximately the maximum packing
#' density that can be achieved in PECC problem.
#'
#' The sample diameter has
#' \eqn{\mu = sample.diameter} and \eqn{\sigma = 0.2} to take into account
#' variations in sample disc preparation (i.e. applying silicon spray to the
#' disc). A lower truncation point at `x = 0.5` is used, which assumes
#' that aliquots with smaller sample diameters of 0.5 mm are discarded.
#' Likewise, the normal distribution is truncated at 9.8 mm, which is the
#' diameter of the sample disc.
#'
#' For each random sample drawn from the
#' normal distributions the amount of grains on the aliquot is calculated. By
#' default, `10^5` iterations are used, but can be reduced/increased with
#' `MC.iter` (see `...`). The results are visualised in a bar- and
#' boxplot together with a statistical summary.
#'
#' @param grain.size [numeric] (**required**):
#' mean grain size (microns) or a range of grain sizes from which the
#' mean grain size is computed (e.g. `c(100,200)`).
#'
#' @param sample.diameter [numeric] (**required**):
#' diameter (mm) of the targeted area on the sample carrier.
#'
#' @param packing.density [numeric] (*with default*):
#' empirical value for mean packing density. \cr
#' If `packing.density = "Inf"` a hexagonal structure on an infinite plane with
#' a packing density of \eqn{0.906\ldots} is assumed.
#'
#' @param MC [logical] (*optional*):
#' if `TRUE` the function performs a Monte Carlo simulation for estimating the
#' amount of grains on the sample carrier and assumes random errors in grain
#' size distribution and packing density. Requires a vector with min and max
#' grain size for `grain.size`. For more information see details.
#'
#' @param grains.counted [numeric] (*optional*):
#' grains counted on a sample carrier. If a non-zero positive integer is provided this function
#' will calculate the packing density of the aliquot. If more than one value is
#' provided the mean packing density and its standard deviation is calculated.
#' Note that this overrides `packing.density`.
#'
#' @param plot [logical] (*with default*):
#' plot output (`TRUE`/`FALSE`)
#'
#' @param ... further arguments to pass (`main, xlab, MC.iter`).
#'
#' @return
#' Returns a terminal output. In addition an
#' [RLum.Results-class] object is returned containing the
#' following element:
#'
#' \item{.$summary}{[data.frame] summary of all relevant calculation results.}
#' \item{.$args}{[list] used arguments}
#' \item{.$call}{[call] the function call}
#' \item{.$MC}{[list] results of the Monte Carlo simulation}
#'
#' The output should be accessed using the function [get_RLum].
#'
#' @section Function version: 0.31
#'
#' @author Christoph Burow, University of Cologne (Germany)
#'
#' @references
#' Duller, G.A.T., 2008. Single-grain optical dating of Quaternary
#' sediments: why aliquot size matters in luminescence dating. Boreas 37,
#' 589-612.
#'
#' Heer, A.J., Adamiec, G., Moska, P., 2012. How many grains
#' are there on a single aliquot?. Ancient TL 30, 9-16.
#'
#' **Further reading**
#'
#' Chang, H.-C., Wang, L.-C., 2010. A simple proof of Thue's
#' Theorem on Circle Packing. [https://arxiv.org/pdf/1009.4322v1](),
#' 2013-09-13.
#'
#' Graham, R.L., Lubachevsky, B.D., Nurmela, K.J.,
#' Oestergard, P.R.J., 1998. Dense packings of congruent circles in a circle.
#' Discrete Mathematics 181, 139-154.
#'
#' Huang, W., Ye, T., 2011. Global
#' optimization method for finding dense packings of equal circles in a circle.
#' European Journal of Operational Research 210, 474-481.
#'
#' @examples
#'
#' ## Estimate the amount of grains on a small aliquot
#' calc_AliquotSize(grain.size = c(100,150), sample.diameter = 1, MC.iter = 100)
#'
#' ## Calculate the mean packing density of large aliquots
#' calc_AliquotSize(grain.size = c(100,200), sample.diameter = 8,
#' grains.counted = c(2525,2312,2880), MC.iter = 100)
#'
#' @md
#' @export
calc_AliquotSize <- function(
grain.size,
sample.diameter,
packing.density = 0.65,
MC = TRUE,
grains.counted,
plot=TRUE,
...
) {
.set_function_name("calc_AliquotSize")
on.exit(.unset_function_name(), add = TRUE)
##==========================================================================##
## CONSISTENCY CHECK OF INPUT DATA
##==========================================================================##
if (missing(grain.size) ||
length(grain.size) == 0 || length(grain.size) > 2) {
.throw_error("Please provide the mean grain size or a range ",
"of grain sizes (in microns)")
}
if(packing.density < 0 | packing.density > 1) {
if(packing.density == "inf") {
} else {
.throw_error("'packing.density' expects values between 0 and 1")
}
}
.validate_positive_scalar(sample.diameter)
if (sample.diameter > 9.8)
.throw_warning("A sample diameter of ", sample.diameter, " mm was ",
"specified, but common sample discs are 9.8 mm in diameter")
if(missing(grains.counted) == FALSE) {
if(MC == TRUE) {
MC = FALSE
cat(paste("\nMonte Carlo simulation is only available for estimating the",
"amount of grains on the sample disc. Automatically set to",
"FALSE.\n"))
}
}
if(MC == TRUE && length(grain.size) != 2) {
.throw_error("'grain.size' must be a vector containing the min and max ",
"grain size when using Monte Carlo simulations")
}
##==========================================================================##
## ... ARGUMENTS
##==========================================================================##
# set default parameters
settings <- list(MC.iter = 10^4,
verbose = TRUE)
# override settings with user arguments
settings <- modifyList(settings, list(...))
##==========================================================================##
## CALCULATIONS
##==========================================================================##
# calculate the mean grain size
range.flag<- FALSE
if(length(grain.size) == 2) {
gs.range<- grain.size
grain.size<- mean(grain.size)
range.flag<- TRUE
}
# use ~0.907... from Thue's Theorem as packing density
if(packing.density == "inf") {
packing.density = pi/sqrt(12)
}
# function to calculate the amount of grains
calc_n<- function(sd, gs, d) {
n<- ((pi*(sd/2)^2)/
(pi*(gs/2000)^2))*d
return(n)
}
# calculate the amount of grains on the aliquot
if(missing(grains.counted) == TRUE) {
n.grains<- calc_n(sample.diameter, grain.size, packing.density)
##========================================================================##
## MONTE CARLO SIMULATION
if(MC == TRUE && range.flag == TRUE) {
# create a random set of packing densities assuming a normal
# distribution with the empirically determined standard deviation of
# 0.18.
d.mc<- rnorm(settings$MC.iter, packing.density, 0.18)
# in a PECC the packing density can not be larger than ~0.87
d.mc[which(d.mc > 0.87)]<- 0.87
d.mc[which(d.mc < 0.25)]<- 0.25
# create a random set of sample diameters assuming a normal
# distribution with an assumed standard deviation of
# 0.2. For a more conservative estimate this is divided by 2.
sd.mc<- rnorm(settings$MC.iter, sample.diameter, 0.2)
# it is assumed that sample diameters < 0.5 mm either do not
# occur, or are discarded. Either way, any smaller sample
# diameter is capped at 0.5.
# Also, the sample diameter can not be larger than the sample
# disc, i.e. 9.8 mm.
sd.mc[which(sd.mc <0.5)]<- 0.5
if (sample.diameter <= 9.8)
sd.mc[which(sd.mc >9.8)]<- 9.8
# create random samples assuming a normal distribution
# with the mean grain size as mean and half the range (min:max)
# as standard deviation. For a more conservative estimate this
# is further devided by 2, so half the range is regarded as
# two sigma.
gs.mc<- rnorm(settings$MC.iter, grain.size, diff(gs.range)/4)
# draw random samples from the grain size spectrum (gs.mc) and calculate
# the mean for each sample. This gives an approximation of the variation
# in mean grain size on the sample disc
gs.mc.sampleMean<- vector(mode = "numeric")
for(i in 1:length(gs.mc)) {
gs.mc.sampleMean[i]<- mean(sample(gs.mc, calc_n(
sample(sd.mc, size = 1),
grain.size,
sample(d.mc, size = 1)
), replace = TRUE))
}
# create empty vector for MC estimates of n
MC.n<- vector(mode="numeric")
# calculate n for each MC data set
for(i in 1:length(gs.mc)) {
MC.n[i]<- calc_n(sd.mc[i],
gs.mc.sampleMean[i],
d.mc[i])
}
# summarize MC estimates
MC.q<- quantile(MC.n, c(0.05,0.95))
MC.n.kde<- density(MC.n, n = 10000)
# apply student's t-test
MC.t.test<- t.test(MC.n)
MC.t.lower<- MC.t.test["conf.int"]$conf.int[1]
MC.t.upper<- MC.t.test["conf.int"]$conf.int[2]
MC.t.se<- (MC.t.upper-MC.t.lower)/3.92
# get unweighted statistics from calc_Statistics() function
MC.stats<- calc_Statistics(as.data.frame(cbind(MC.n,0.0001)))$unweighted
}
}#EndOf:estimate number of grains
##========================================================================##
## CALCULATE PACKING DENSITY
if(missing(grains.counted) == FALSE) {
area.container<- pi*sample.diameter^2
if(length(grains.counted) == 1) {
area.grains<- (pi*(grain.size/1000)^2)*grains.counted
packing.density<- area.grains/area.container
}
else {
packing.densities<- length(grains.counted)
for(i in 1:length(grains.counted)) {
area.grains<- (pi*(grain.size/1000)^2)*grains.counted[i]
packing.densities[i]<- area.grains/area.container
}
std.d<- sd(packing.densities)
}
}
##==========================================================================##
##TERMINAL OUTPUT
##==========================================================================##
if (settings$verbose) {
cat("\n [calc_AliquotSize]")
cat(paste("\n\n ---------------------------------------------------------"))
cat(paste("\n mean grain size (microns) :", grain.size))
cat(paste("\n sample diameter (mm) :", sample.diameter))
if(missing(grains.counted) == FALSE) {
if(length(grains.counted) == 1) {
cat(paste("\n counted grains :", grains.counted))
} else {
cat(paste("\n mean counted grains :", round(mean(grains.counted))))
}
}
if(missing(grains.counted) == TRUE) {
cat(paste("\n packing density :", round(packing.density,3)))
}
if(missing(grains.counted) == FALSE) {
if(length(grains.counted) == 1) {
cat(paste("\n packing density :", round(packing.density,3)))
} else {
cat(paste("\n mean packing density :", round(mean(packing.densities),3)))
cat(paste("\n standard deviation :", round(std.d,3)))
}
}
if(missing(grains.counted) == TRUE) {
cat(paste("\n number of grains :", round(n.grains,0)))
}
if(MC == TRUE && range.flag == TRUE) {
cat(paste(cat(paste("\n\n --------------- Monte Carlo Estimates -------------------"))))
cat(paste("\n number of iterations (n) :", settings$MC.iter))
cat(paste("\n median :", round(MC.stats$median)))
cat(paste("\n mean :", round(MC.stats$mean)))
cat(paste("\n standard deviation (mean) :", round(MC.stats$sd.abs)))
cat(paste("\n standard error (mean) :", round(MC.stats$se.abs, 1)))
cat(paste("\n 95% CI from t-test (mean) :", round(MC.t.lower), "-", round(MC.t.upper)))
cat(paste("\n standard error from CI (mean):", round(MC.t.se, 1)))
cat(paste("\n ---------------------------------------------------------\n"))
} else {
cat(paste("\n ---------------------------------------------------------\n"))
}
}
##==========================================================================##
##RETURN VALUES
##==========================================================================##
# prepare return values for mode: estimate grains
if(missing(grains.counted) == TRUE) {
summary<- data.frame(grain.size = grain.size,
sample.diameter = sample.diameter,
packing.density = packing.density,
n.grains = round(n.grains,0),
grains.counted = NA)
}
# prepare return values for mode: estimate packing density/densities
if(missing(grains.counted) == FALSE) {
# return values if only one value for counted.grains is provided
if(length(grains.counted) == 1) {
summary<- data.frame(grain.size = grain.size,
sample.diameter = sample.diameter,
packing.density = packing.density,
n.grains = NA,
grains.counted = grains.counted)
} else {
# return values if more than one value for counted.grains is provided
summary<- data.frame(rbind(1:5))
colnames(summary)<- c("grain.size", "sample.diameter", "packing.density",
"n.grains","grains.counted")
for(i in 1:length(grains.counted)) {
summary[i,]<- c(grain.size, sample.diameter, packing.densities[i],
n.grains = NA, grains.counted[i])
}
}
}
if(!MC) {
MC.n<- NULL
MC.stats<- NULL
MC.n.kde<- NULL
MC.t.test<- NULL
MC.q<- NULL
}
if(missing(grains.counted)) grains.counted<- NA
call<- sys.call()
args<- as.list(sys.call())[-1]
# create S4 object
newRLumResults.calc_AliquotSize <- set_RLum(
class = "RLum.Results",
data = list(
summary=summary,
MC=list(estimates=MC.n,
statistics=MC.stats,
kde=MC.n.kde,
t.test=MC.t.test,
quantile=MC.q)),
info = list(call=call,
args=args))
##=========##
## PLOTTING
if(plot==TRUE) {
try(plot_RLum.Results(newRLumResults.calc_AliquotSize, ...))
}
# Return values
invisible(newRLumResults.calc_AliquotSize)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.