# create and plot data set that consists of five standards with concentrations 
# of 20, 40, 60, 80, and 100 ppb; A = ebC = kC with k = 0.012 and with a random
# uniform offset to provide some scatter in calibration data; data collection is
# 10 seconds of a reagent blank or a standard bounded on either side by 10 sec of
# flushing system with water (mimicing automated collection of data by FIA);
# sampling rate is 10 points per second, or 100 points per aquisition period;
# blank is the middle third of gap0, which is 10 sec flush, 10 sec blank; 10 sec
# flush; random normal noise with mean of 0 and sd of 0.025 included
concstd = seq(20, 100, 20)
k = 0.012
abs = k * concstd + runif(5, -0.025, +0.025)
gap0 = rnorm(300, 0, 0.025)
std1 = abs[1] + rnorm(100, 0, 0.025)
gap1 = rnorm(100, 0, 0.025)
std2 = abs[2] + rnorm(100, 0, 0.025)
gap2 = rnorm(100, 0, 0.025)
std3 = abs[3] + rnorm(100, 0, 0.025)
gap3 = rnorm(100, 0, 0.025)
std4 = abs[4] + rnorm(100, 0, 0.025)
gap4 = rnorm(100, 0, 0.025)
std5 = abs[5] + rnorm(100, 0, 0.025)
gap5 = rnorm(101, 0, 0.025)
signal = c(gap0, std1, gap1, std2, gap2, std3, gap3, std4, gap4, std5, gap5)
time = seq(0, 130, 0.1)
plot(time, signal, type = "l", ylim = c(-0.1, 1.4), xlab = "time (s)")
abline(v = 10.1, lty = 3)
abline(v = 20.0, lty = 3)
text(x = 15, y = 0.2, srt = 90, "blank")
text(x = 35, y = 0.37, "20 ppb")
text(x = 55, y = 0.60, "40 ppb")
text(x = 75, y = 0.85, "60 ppb")
text(x = 95, y = 1.08, "80 ppb")
text(x = 115, y = 1.33, "100 ppb")

# extract means for blank and calibration standards; plot calibration data; build linear model and display slope and intercept
sigstand = c(mean(signal[101:200]), mean(signal[301:400]), mean(signal[501:600]), mean(signal[701:800]), mean(signal[901:1000]), mean(signal[1101:1200]))
conc = c(0, concstd)
plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
model = lm(sigstand ~ conc)
abline(model)
summary(model)
legend(x = "topleft", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 6))), bty = "n")

sd.blank = sd(signal[101:200])
LOD = (3 * sd.blank - model$coefficients[[1]])/model$coefficients[[2]]
LOQ = (10 * sd.blank - model$coefficients[[1]])/model$coefficients[[2]]
plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
model = lm(sigstand ~ conc)
abline(model)
legend(x = "topleft", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 6)), paste0("sd of blank: ", round(sd.blank, digits = 4)), paste0("LOD: ", round(LOD, digits = 2), " ppb"), paste0("LOQ: ", round(LOQ, digits = 2), " ppb")), bty = "n")
lines(x = c(-10, LOD), y = c(3 * sd.blank, 3 * sd.blank), lty = 2)
lines(x = c(-10, LOQ), y = c(10 * sd.blank, 10 * sd.blank), lty = 2)
lines(x = c(LOD, LOD), y = c(-1, 3 * sd.blank), lty = 2)
lines(x = c(LOQ, LOQ), y = c(-1, 10 * sd.blank), lty = 2)
text(x = LOD, y = 0, "LOD", pos = 4, offset = 0.1)
text(x = LOQ, y = 0.1, "LOQ", pos = 4, offset = 0.1)

# add interferent with concentration of 2000 ppb; plot signal
# int = (0.008 * 20 + runif(1, -0.025, +0.025)) + rnorm(100, 0, 0.025)
int = (0.008 * 20)
sigint = c(gap0[1:100], gap0[101:200] + int, gap0[201:300], std1 + int, gap1, std2 + int, gap2, std3 + int, gap3, std4 + int, gap4, std5 + int, gap5)
plot(time, sigint, type = "l",  xlab = "time (s)", ylim = c(0, 1.6))
abline(v = 10.1, lty = 3)
abline(v = 20.0, lty = 3)
text(x = 15, y = 0.40, srt = 90, "blank")
text(x = 35, y = 0.60, "20 ppb")
text(x = 55, y = 0.85, "40 ppb")
text(x = 75, y = 1.05, "60 ppb")
text(x = 95, y = 1.30, "80 ppb")
text(x = 115, y = 1.55, "100 ppb")

# plot calibration curve
sigstand = c(mean(sigint[101:200]), mean(sigint[301:400]), mean(sigint[501:600]), mean(sigint[701:800]), mean(sigint[901:1000]), mean(sigint[1101:1200]))
conc = c(0, concstd)
plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
model = lm(sigstand ~ conc)
abline(model)
sd.blank = sd(sigint[101:200])
legend(x = "topleft", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 4)), paste0("sd of blank: ", round(sd.blank, digits = 3)), paste0("[interferent]: 2000 ppb")), bty = "n")

# calculate k for interferent using intercept, and K_A/I
k.i = model$coefficients[1]/2000
k.a = k
K.ai = k.i/k.a
plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
abline(model)
legend(x = "topleft", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 4)), paste0("sd of blank: ", round(sd.blank, digits = 3)), paste0("[interferent]: 2000 ppb"), paste0("k_analyte: ", round(model$coefficients[2], digits = 4)), paste0("k_interferent: ", round(k.i, digits = 7)), paste0("K_A/I: ", round(K.ai, digits = 6))), bty = "n")
# legend(x = "topleft", legend = c(paste0("k_analyte: ", round(model$coefficients[2], digits = 4)), paste0("k_interferent: ", round(model$coefficients[1], digits = 4)), paste0("K_A/I: ", round(K.ai, digits = 6))), bty = "n")

# calculate LOD and LOQ, report values, and add lines to plot
a.lod = 3 * sd.blank
a.loq = 10 * sd.blank
i.lod = a.lod/k.i
i.loq = a.loq/k.i
plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
abline(model)
legend(x = "bottomright", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 4)), paste0("sd of blank: ", round(sd.blank, digits = 3)), paste0("[interferent]: 2000 ppb"), paste0("k_analyte: ", round(model$coefficients[2], digits = 4)), paste0("k_interferent: ", round(k.i, digits = 7)), paste0("K_A/I: ", round(K.ai, digits = 6)), paste0("[I] equal to LOD: ", round(i.lod, digits = 0), " ppb"), paste0("[I] equal to LOQ: ", round(i.loq, digits = 0), " ppb")), bty = "n")

# LOD = (3 * sd.blank - model$coefficients[[1]])/model$coefficients[[2]]
# LOQ = (10 * sd.blank - model$coefficients[[1]])/model$coefficients[[2]]
# plot(conc, sigstand, pch = 19, ylim = c(0, 1.4), xlim = c(0, 100), xlab = "concentration (ppb)", ylab = "average signal")
# model = lm(sigstand ~ conc)
# abline(model)
# legend(x = "topleft", legend = c(paste0("slope: ", round(model$coefficients[2], digits = 4)), paste0("intercept: ", round(model$coefficients[1], digits = 4)), paste0("sd of blank: ", round(sd.blank, digits = 4)), paste0("LOD: ", round(LOD, digits = 2), " ppb"), paste0("LOQ: ", round(LOQ, digits = 2), " ppb")), bty = "n")
# lines(x = c(-10, LOD), y = c(3 * sd.blank, 3 * sd.blank), lty = 2)
# lines(x = c(-10, LOQ), y = c(10 * sd.blank, 10 * sd.blank), lty = 2)
# lines(x = c(LOD, LOD), y = c(-1, 3 * sd.blank), lty = 2)
# lines(x = c(LOQ, LOQ), y = c(-1, 10 * sd.blank), lty = 2)
# text(x = LOD, y = 0, "LOD", pos = 4, offset = 0.1)
# text(x = LOQ, y = 0.1, "LOQ", pos = 4, offset = 0.1)