Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
RLumSK opened this issue Jul 3, 2020 · 0 comments
Open

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

RLumSK opened this issue Jul 3, 2020 · 0 comments
Assignees

Comments

@RLumSK
Copy link
Member

RLumSK commented Jul 3, 2020

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)?

@RLumSK RLumSK added the bug label Jul 3, 2020
@RLumSK RLumSK self-assigned this Jul 3, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

1 participant