# application of beer's law

allSpec <- read_csv("allSpec.csv")
load("co_stds.RData")
load("cr_stds.RData")
load("cu_stds.RData")
load("ni_stds.RData")

# external standardization for copper

# examine spectra for cu stds
matplot(cu_stds2$wavelength, cu_stds2[2:6], type = "l", lwd = 2, xlab = "wavelength (nm)", ylab = "absorbance")
legend("topleft", legend = c("0.0500 M", "0.0400 M", "0.0300 M", "0.0200 M", "0.0100 M"), lwd = 2, lty = 1:5, col = 1:5, bty = "n")

# find the wavelength where absorbance is a maximum
cu.id = which.max(cu_stds2$cu_std1)
abline(v = cu_stds2$wavelength[cu.id], lwd = 2, lty = 2)
lambda = cu_stds2$wavelength[cu.id]
lambda

# find absorbances at this wavelength
cu_abs = cu_stds2[cu.id, 2:6]
cu_abs

# plot the calibration curve
plot(x = cu_stdconc, y = cu_abs, pch = 19, col = "blue", xlab = "[Cu] (M)", ylab = "absorbance", xlim = c(0, max(cu_stdconc)), ylim = c(0, max(cu_stds2[cu.id,2])))
grid()

# construct the linear model for the calibration curve
cu.lm = lm(as.numeric(cu_stds2[cu.id, 2:6]) ~ cu_stdconc)
summary(cu.lm)
abline(cu.lm, lwd = 2, col = "blue", lty = 2)

# determine concentration of copper in unknown
cu_unkabs = rnorm(3, 0.3, 0.01)
cu_unkabs
library(chemCal)
cu_unkconc = inverse.predict(cu.lm, cu_unkabs)
cu_unkconc

# external standardization for copper and nickel

# examine the spectra for copper and nickel
wavelengths = as.numeric(colnames(allSpec[8:642]))
plot(x = cu_stds2$wavelength, y = cu_stds2$cu_std1, xlab = "wavelength (nm)", ylab = "absorbance", type = "l", lwd = 2, lty = 2, col = "blue")
lines(x = ni_stds2$wavelength, y = ni_stds2$ni_std1, lty = 2, lwd = 2, col = "forestgreen")
legend("top", legend = c("", "Cu standard", "Ni standard"), lwd = 2, lty = c(1, 2, 2), col = c("white", "blue", "forestgreen"), bty = "n")
grid()

# find wavelength of maximum absorbance for copper and for nickel
cu.id = which.max(cu_stds2$cu_std1)
ni.id = which.max(ni_stds2$ni_std1)
lambda_cu = cu_stds2$wavelength[cu.id]
lambda_ni = ni_stds2$wavelength[ni.id]
lambda_cu
lambda_ni
abline(v = cu_stds2$wavelength[cu.id], lwd = 2, lty = 3, col = "blue")
abline(v = ni_stds2$wavelength[ni.id], lwd = 2, lty = 3, col = "forestgreen")

# examine spectrum for mixture
lines(x = wavelengths, y = allSpec[62, 8:642], lty = 1, lwd = 2, col = "black")
legend("top", legend = c("Cu/Ni Mixture", "", ""), lwd = 2, lty = c(1, 2, 2), col = c("black", "blue", "forestgreen"), bty = "n")

# plot calibration curves for copper and nickel
plot(x = cu_stdconc, y = cu_stds2[cu.id,2:6], pch = 19, col = "blue", xlab = "[Cu] or [Ni] (M)", ylab = "absorbance", xlim = c(0, max(ni_stdconc)), ylim = c(0, max(cu_stds2[cu.id,2])))
cu.lm1 = lm(as.numeric(cu_stds2[cu.id, 2:6]) ~ cu_stdconc)
abline(cu.lm1, lwd = 2, col = "blue", lty = 2)
points(x = cu_stdconc, y = cu_stds2[ni.id,2:6], pch = 15, col = "blue")
cu.lm2 = lm(as.numeric(cu_stds2[ni.id, 2:6]) ~ cu_stdconc)
abline(cu.lm2, lwd = 2, col = "blue", lty = 2)
points(x = ni_stdconc, y = ni_stds2[cu.id,2:6], pch = 19, col = "forestgreen")
ni.lm1 = lm(as.numeric(ni_stds2[cu.id, 2:6]) ~ ni_stdconc)
abline(ni.lm1, lwd = 2, col = "forestgreen", lty = 2)
points(x = ni_stdconc, y = ni_stds2[ni.id,2:6], pch = 15, col = "forestgreen")
ni.lm2 = lm(as.numeric(ni_stds2[ni.id, 2:6]) ~ ni_stdconc)
abline(ni.lm2, lwd = 2, col = "forestgreen", lty = 2)
legend("right", legend = c("Cu: 809.1 nm", "Cu: 394.2 nm", "Ni: 809.1 nm", "Ni: 394.2 nm"), pch = c(19,15,19,15), col = c("blue", "blue", "forestgreen", "forestgreen"), bty = "n")

# obtain eb values for copper and nickel at each wavelength as a matrix
eb = matrix(data = c(eb_cu1 = cu.lm1$coefficients[2], 
                     eb_cu2 = cu.lm2$coefficients[2], 
                     eb_ni1 =  ni.lm1$coefficients[2], 
                     eb_ni2 = ni.lm2$coefficients[2]), 
            nrow = 2, ncol = 2, byrow = FALSE)
colnames(eb) = c("Cu", "Ni")
rownames(eb) = c("809.1 nm", "394.2 nm")
round(eb, digits = 3)

# obtain absorbance values for copper and nickel at each wavelength
abs = matrix(data = c(allSpec[62, cu.id], allSpec[62, ni.id]), 
             nrow = 2, ncol = 1, byrow = FALSE)
colnames(abs) = "mixture"
rownames(abs) = c("809.1 nm", "394.2 nm")
round(as.numeric(abs), digits = 3)

# solve the simultaneous equations
mix_conc = solve(eb, abs)
mix_conc

# actual concentrations
cu_actual = allSpec[62, "concCu"]
cu_actual
cu_error = 100 * (cu_actual - mix_conc[1])/cu_actual
round(cu_error, digits = 2)

ni_actual = allSpec[62, "concNi"]
ni_actual
ni_error = 100 * (ni_actual - mix_conc[2])/ni_actual
round(ni_error, digits = 2)