old.par = par(mfrow = c(2, 2))
# experimental retention times as a function of pH
tr.ph = c(13.36, 11.57, 10.08, 8.43, 7.04, 5.60, 4.08)
tr.ba = c(39.68, 29.74, 23.58, 18.41, 13.88, 11.01, 7.65)
tr.tp = c(16.10, 7.73, 4.45, 2.92, 1.84, 1.33, 0.93)
tr.pa = c(9.10, 8.54, 8.07, 7.12, 6.45, 5.42, 4.11)
ph = c(3.65, 4.00, 4.25, 4.42, 4.63, 4.82, 5.01)
plot(ph, tr.ba, pch = 19, col = "blue", ylim = c(0, max(tr.ba)), xlim = c(3.5, 5.5), xlab = "pH", ylab = "retention time")
points(ph, tr.pa, pch = 19, col = "red")
points(ph, tr.ph, pch = 19, col = "green")
points(ph, tr.tp, pch = 19, col = "orange")
legend(x = "topright", legend = c("BA", "PA", "PH", "TP"), col = c("blue", "red", "green", "orange"), lty = 1, bty = "n")

# predicted retention times determined by best fit
pH = seq(3.5, 5.5, 0.01)
r = 10^(pH - 4.71)
pred.ba = (53.07 + 2.76*4.25*r)/(1+4.25*r)
pred.pa = (9.61 + 0.67*0.68*r)/(1+0.68*r)
pred.ph = (15.02 + 0.75*1.58*r)/(1+1.58*r)
pred.tp = (156 + 0.076*100.1*r)/(1+100.1*r)
plot(ph, tr.ba, pch = 19, col = "blue", ylim = c(0, max(tr.ba)), xlim = c(3.5, 5.5), xlab = "pH", ylab = "retention time")
points(ph, tr.pa, pch = 19, col = "red")
points(ph, tr.ph, pch = 19, col = "green")
points(ph, tr.tp, pch = 19, col = "orange")
legend(x = "topright", legend = c("BA", "PA", "PH", "TP"), col = c("blue", "red", "green", "orange"), lty = 1, bty = "n")
lines(pH, pred.ba, col = "blue")
lines(pH, pred.pa, col = "red")
lines(pH, pred.ph, col = "green")
lines(pH, pred.tp, col = "orange")

# predicted separation factors from predicted retention times
ba.pa = abs((pred.ba - pred.pa))/(pred.ba + pred.pa)
ba.ph = abs((pred.ba - pred.ph))/(pred.ba + pred.ph)
ba.tp = abs((pred.ba - pred.tp))/(pred.ba + pred.tp)
pa.ph = abs((pred.pa - pred.ph))/(pred.pa + pred.ph)
pa.tp = abs((pred.pa - pred.tp))/(pred.pa + pred.tp)
ph.tp = abs((pred.ph - pred.tp))/(pred.ph + pred.tp)
plot(pH, ba.pa, type = "l", lty = 1, lwd = 2, ylab = "separation factor", ylim = c(0, 0.25), xlab = "pH")
lines(pH, ba.ph, lty = 2, lwd = 2)
lines(pH, ba.tp, lty = 3, lwd = 2)
lines(pH, pa.ph, lty = 4, lwd = 2)
lines(pH, pa.tp, lty = 5, lwd = 2)
lines(pH, ph.tp, lty = 6, lwd = 2)
legend(x = "topright", legend = c("BA & PA", "BA & PH", "BA & TP", "PA & PH", "PA & TP", "PH & TP"), lty = c(1, 2, 3, 4, 5, 6), lwd = 2, bty = "n")

# plot window diagram
plot(pH, ba.pa, type = "l", lty = 1, lwd = 2, ylab = "separation factor", ylim = c(0, 0.25), xlab = "pH")
lines(pH, ba.ph, lty = 2, lwd = 2)
lines(pH, ba.tp, lty = 3, lwd = 2)
lines(pH, pa.ph, lty = 4, lwd = 2)
lines(pH, pa.tp, lty = 5, lwd = 2)
lines(pH, ph.tp, lty = 6, lwd = 2)
legend(x = "topright", legend = c("BA & PA", "BA & PH", "BA & TP", "PA & PH", "PA & TP", "PH & TP"), lty = c(1, 2, 3, 4, 5, 6), lwd = 2, bty = "n")
window = pmin(ba.pa, ba.ph, ba.tp, pa.ph, pa.tp, ph.tp)
dx = seq(3.5, 5.5, 0.01)
polygon(c(3.5, dx, 5.5), y = c(0, window, 0), border = NA, col = "gray")
par(old.par)