R-Lum/Luminescence

calc_Lamothe2003() - wrong fading corrected Ln/Tn error?

Opened this issue · 1 comments

Transferred issue from an e-mail discussion with the Lux group in Montréal, for the example below the Ln/Tn error after fading correction appears to be too small.

library(Luminescence)
df <- data.frame(
  Dose = c(0,600,150,300,1200,3200,0,600),
  LxTx = c(2.124,2.29,0.613,1.166,3.879,7.857,0.021,2.184),
  LxTx_X = c(0.046,0.049,0.013,0.025,0.083,0.169,0.001,0.047)
)

results <- calc_Lamothe2003(
  object = df,
  dose_rate.envir = c(1.75,0.02),
  dose_rate.source = c(0.125,0.002),
  g_value = c(4.0,0.2),
  tc = 400,
  tc.g_value = 172800,
  plot = FALSE
)

cat(">> Original LnTn error:", round(df$LxTx_X[1]/df$LxTx[1] * 100,1), "%","\n")
cat(">> Corrected LnTn error:", round(results$data$LnTn_AFTER.ERROR/results$data$LnTn_AFTER * 100,1), "%")

The error calculation might be uncorrected, in particular the following part needs to be checked:

## fading correction (including dose rate conversion from Gy/s to Gy/ka)
## and error calculation
## the formula in Lamothe et al. (2003) reads:
## I_faded = I_unfaded*(1-g*log((1/e)*DR_lab/DR_soil)))
rr <- 31.5576e+09 * dose_rate.source[1] / (exp(1) * dose_rate.envir[1])
s_rr <- sqrt((dose_rate.source[2]/dose_rate.source[1])^2 + (dose_rate.envir[2]/dose_rate.envir[1])^2) * rr
Fading_C <- 1 - g_value[1] / 100 * log10(rr)
sFading_C <- sqrt((log10(rr) * g_value[2]/100)^2 + (g_value[1]/(100 * rr) * s_rr)^2)
# store original Lx/Tx in new object
LnTn_BEFORE <- data[[2]][1]
LnTn_BEFORE.ERROR <- data[[3]][1]
# apply to input data
data[[2]][1] <- data[[2]][1] / Fading_C
data[[3]][1] <- sqrt((data[[3]][1]/data[[2]][1])^2 +
((1/Fading_C - 1) * sFading_C/Fading_C)^2) * data[[2]][1]

Why we have the term (1/Fading_C - 1)?

@mcol Can you double check the math here?