# load data
load("~/Box Sync/p-harvey/Teaching/Chem 450/Class Materials/stdMethods.RData")

# plot spectra for five external standards, a calibration std, and a sample
plot(wavelength, extstd5, type = "l", lwd = 2, xlab = "wavelength", ylab = "absorbance")
lines(wavelength, extstd4, lwd = 2)
lines(wavelength, extstd3, lwd = 2)
lines(wavelength, extstd2, lwd = 2)
lines(wavelength, extstd1, lwd = 2)
lines(wavelength, calstd, lwd = 2, lty = 2)
lines(wavelength, sample1, lwd = 2, lty = 3)
legend(x = "topright", legend = c("external standard", "calibration standard", "sample"), lty = c(1, 2, 3), bty = "n")

# identify analytical wavelength
id = which.max(extstd5)
wavelength[id]
abline(v = wavelength[id], lwd = 2, lty = 2)

# build linear regression model for external standards
A.extstd = c(extstd1[id], extstd2[id], extstd3[id], extstd4[id], extstd5[id])
C.extstd = c(0.20, 0.40, 0.60, 0.80, 1.0)
extstd.lm = lm(A.extstd ~ C.extstd)
summary(extstd.lm)
plot(C.extstd, A.extstd, pch = 19, xlab = "concentration of external standards (M)", ylab = "absorbance", ylim = c(0, 0.35), xlim = c(0, 1))
abline(extstd.lm)
legend(x = "topleft", legend = c(paste0("slope: ", round(extstd.lm$coefficients[2], digits = 4)), paste0("intercept: ", round(extstd.lm$coefficients[1], digits = 4))), bty = "n")

# evaluate calibration model using calibration standard with C = 0.50 M
A.calstd = calstd[id]
A.calstd
abline(h = A.calstd, lty = 4)
text(x = 0.8, y = 0.17, labels = paste0("A = ", round(A.calstd, digits = 3)))

C.calstd = (A.calstd - extstd.lm$coefficients[[1]])/extstd.lm$coefficients[[2]]
abline(v = C.calstd, lty = 4)
round(C.calstd, digits = 3)
per.error = 100 * (C.calstd - 0.50)/0.50
noquote(paste0(round(per.error, digits = 2), "% error"))

# evaluate sample with C = ???
A.sample1 = sample1[id]
A.sample1
abline(h = A.sample1, lty = 5)
text(x = 0.8, y = 0.13, labels = paste0("A = ", round(A.sample1, digits = 3)))

C.sample1 = (A.sample1 - extstd.lm$coefficients[[1]])/extstd.lm$coefficients[[2]]
abline(v = C.sample1, lty = 5)
round(C.sample1, digits = 3)

# evaluate sample using spike recovery
# 5 mL of sample + 1 mL of 1.0 M ext std diluted to 10 mL
A.spike = spike[id]
A.spike
abline(h = A.spike, lty = 6)
text(x = 0.8, y = 0.095, labels = paste0("A = ", round(A.spike, digits = 3)))

C.spike = (A.spike - extstd.lm$coefficients[[1]])/extstd.lm$coefficients[[2]]
abline(v = C.spike, lty = 6)
round(C.spike, digits = 3)

C.added = 1.0 * 0.1
C.added

C.samdiluted = 0.5 * C.sample1
round(C.samdiluted, digits = 3)

per.recov = 100 * (C.spike - C.samdiluted)/C.added
noquote(paste0(round(per.recov, digits = 1), "% error"))

# standard additions: 0, 1, 2, 3, 4 mL of 1.0 M and 5 mL sample diluted to 10 mL
plot(wavelength, stdadd4, type = "l", lwd = 2, xlab = "wavelength", ylab = "absorbance", lty = 2)
lines(wavelength, stdadd3, lwd = 2, lty = 2)
lines(wavelength, stdadd2, lwd = 2, lty = 2)
lines(wavelength, stdadd1, lwd = 2, lty = 2)
lines(wavelength, stdadd0, lwd = 2, lty = 1)
abline(v = wavelength[id], lwd = 2, lty = 2)
legend(x = "topright", legend = c("original sample (diluted 5:10)", "with standard additions"), lty = c(1, 2), bty = "n")

A.stdadd = c(stdadd0[id], stdadd1[id], stdadd2[id], stdadd3[id], stdadd4[id])
C.stdadd = c(0, 0.10, 0.20, 0.30, 0.40)
stdadd.lm = lm(A.stdadd ~ C.stdadd)
summary(stdadd.lm)
plot(C.stdadd, A.stdadd, pch = 19, xlab = "concentration added (M)", ylab = "absorbance", ylim = c(0, 0.2), xlim = c(-0.5, 0.5))
abline(v = 0)
abline(stdadd.lm)
legend(x = "topleft", legend = c(paste0("slope: ", round(stdadd.lm$coefficients[2], digits = 4)), paste0("intercept: ", round(stdadd.lm$coefficients[1], digits = 4))), bty = "n")

C.sample = 2 * stdadd.lm$coefficients[[1]]/stdadd.lm$coefficients[[2]]
round(C.sample, digits = 3)
per.error = 100 * (C.sample - 0.50)/0.50
noquote(paste0(round(per.error, digits = 2), "% error"))

# what if we cannot control pathlength?
plot(wavelength, internal5, type = "l", lwd = 2, xlab = "wavelength", ylab = "absorbance")
lines(wavelength, internal4, lwd = 2)
lines(wavelength, internal3, lwd = 2)
lines(wavelength, internal2, lwd = 2)
lines(wavelength, internal1, lwd = 2)
abline(v = wavelength[id], lwd = 2, lty = 2)

A.internal = c(internal1[id], internal2[id], internal3[id], internal4[id], internal5[id])
C.internal = c(0.20, 0.40, 0.60, 0.80, 1.0)
internal.lm = lm(A.internal ~ C.internal)
summary(internal.lm)
plot(C.internal, A.internal, pch = 19, xlab = "concentration (M)", ylab = "absorbance", ylim = c(0, 0.35), xlim = c(0, 1))
abline(internal.lm)

# internal standards: 5 mL of ext std + 1 mL 1.0 M internal std diluted to 10 mL
plot(wavelength, intstd5, type = "l", lwd = 2, xlab = "wavelength", ylab = "absorbance")
lines(wavelength, intstd4, lwd = 2)
lines(wavelength, intstd3, lwd = 2)
lines(wavelength, intstd2, lwd = 2)
lines(wavelength, intstd1, lwd = 2)
abline(v = wavelength[id], lwd = 2, lty = 2)
abline(v = wavelength[521], lwd = 2, lty = 2)

A.internal = c(intstd1[id]/intstd1[521], intstd2[id]/intstd2[521], intstd3[id]/intstd3[521], intstd4[id]/intstd4[521], intstd5[id]/intstd5[521])
C.internal = c(0.20, 0.40, 0.60, 0.80, 1.0)
internal.lm = lm(A.internal ~ C.internal)
summary(internal.lm)
plot(C.internal, A.internal, pch = 19, xlab = "concentration (M)", ylab = "Asample/Aintstd", ylim = c(0, 3), xlim = c(0, 1))
abline(internal.lm)
legend(x = "topleft", legend = c(paste0("slope: ", round(internal.lm$coefficients[2], digits = 3)), paste0("intercept: ", round(internal.lm$coefficients[1], digits = 3))), bty = "n")

S.sample2 = sample2[id]/sample2[521]
S.sample2
abline(h = S.sample2, lty = 4)
text(x = 0.8, y = 1.3, labels = paste0("A = ", round(S.sample2, digits = 3)))


C.sample2 = (S.sample2 - internal.lm$coefficients[[1]])/internal.lm$coefficients[[2]]
abline(v = C.sample2, lty = 4)
round(C.sample2, digits = 3)
per.error = 100 * (C.sample2 - 0.50)/0.50
noquote(paste0(round(per.error, digits = 2), "% error"))