#' Enter the number of extractions to explore, and the total volume for the aqueous phase and the organic phase. For the organic phase, the total volume is divided evenly among the individual extractions.
n = 1
v.aq = 50
v.org = 50/n

#' Enter the pKa values, the Kd values, the initial concentrations, and the selectivity coefficient for the analyte and the interferent

pka.a = 3
pka.i = 8
kd.a = 50
kd.i = 100
conc.a = 1
conc.i = 10
k.ai = 0.5

#' Set up range of pH values to explore and convert to concentration of H3O+

ph = seq(1, 14, 0.01)
h = 10^-ph

#' Calculate distribution ratios assuming analyte is weak acid and interferent is weak base

ka.a = 10^-pka.a
ka.i = 10^-pka.i
d.a = kd.a * h/(h + ka.a)
d.i = kd.i * ka.i/(h + ka.i)

#' Calculate the fraction of analyte and interferent in each phase; these values are equivalent to recoveries.
qaq.a = (v.aq/(d.a*v.org + v.aq))^n
qaq.i = (v.aq/(d.i*v.org + v.aq))^n
qorg.a = 1 - qaq.a
qorg.i = 1 - qaq.i

#' Calculate and report the error for the recovery of the analyte in each phase.
erroraq = 100 * ((qaq.a - 1) + k.ai * conc.i * qaq.i/conc.a)
errororg = 100 * ((qorg.a - 1) + k.ai * conc.i * qorg.i/conc.a)

#' Examine results as a function of pH
plot(ph, qaq.a, type = "l", lwd = 2, lty = 1, col = "blue", 
     xlab = "pH", ylim = c(0, 1), ylab = "fraction")
lines(ph, qaq.i, type = "l", lwd = 2, lty = 1, col = "red")
lines(ph, qorg.a, lwd = 2, lty = 2, col = "blue")
lines(ph, qorg.i, lwd = 2, lty = 2, col = "red")
legend(x = "right", 
       legend = c("Analyte(aq)", "Interferent(aq)", "Analyte(org)", "Interferent(org)"), 
       lwd = 2, lty = c(1, 1, 2, 2), col = c("blue", "red", "blue", "red"), bty = "n")
grid(col = "black")

plot(ph, erroraq, type = "l", lwd = 2, lty = 1, col = "blue", 
     xlab = "pH", ylab = "percent error", ylim = c(-100, 500))
lines(ph, errororg, lwd = 2, lty = 2, col = "blue")
legend(x = "right", 
       legend = c("Analyte(aq)", "Analyte(org)"), lwd = 2, 
       lty = c(1, 2), col = "blue", bty = "n")
grid(col = "black")