# 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 FAAS);
# 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 the blank and the 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 = 4))), bty = "n")

# calculate and display sd for blank
sd.blank = sd(signal[101:200])
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))), bty = "n")

# calculate LOD and LOQ, report values, and add lines to plot
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)