ccextract = function(pH){
  steps = 50
  h = 10^-pH
  kd.a = 5
  kd.b = 4
  kd.c = 6
  ka1.a = 1e-5
  ka2.a = 1e-9
  ka1.b = 1e-3
  ka2.b = 1e-7
  ka1.c = 1e-7
  ka2.c = 1e-11
  f.a = ka1.a * h/(h^2 + ka1.a * h + ka1.a * ka2.a)
  f.b = ka1.b * h/(h^2 + ka1.b * h + ka1.b * ka2.b)
  f.c = ka1.c * h/(h^2 + ka1.c * h + ka1.c * ka2.c)
  d.a = kd.a * f.a
  d.b = kd.b * f.b
  d.c = kd.c * f.c
  q.a = 1/(d.a + 1)
  q.b = 1/(d.b + 1)
  q.c = 1/(d.c + 1)
  p.a = 1 - q.a
  p.b = 1 - q.b
  p.c = 1 - q.c
  f.a = matrix(0, steps, steps)
  f.b = matrix(0, steps, steps)
  f.c = matrix(0, steps, steps)
  for (n in 1:steps){
    for (r in 1:n){
      f.a[r, n] = 
        (factorial(n)/(factorial(n-r)*factorial(r)))*p.a^r*q.a^(n - r)
      f.b[r, n] = 
        (factorial(n)/(factorial(n-r)*factorial(r)))*p.b^r*q.b^(n - r)
      f.c[r, n] = 
        (factorial(n)/(factorial(n-r)*factorial(r)))*p.c^r*q.c^(n - r)
    }
  }
  out = list("fractionA" = f.a, "fractionB" = f.b, "fractionC" = f.c, 
             "steps" = steps, "ph" = pH)
}


cce.plot = function(x){
  plot(x$fractionA[, x$steps], 
       ylim = c(0, max(c(x$fractionA[ , x$steps], 
                         x$fractionB[, x$steps],
                         x$fractionC[, x$steps],
                         1 - sum(x$fractionA[, x$steps]), 
                         1 - sum(x$fractionB[, x$steps]),
                         1 - sum(x$fractionC[, x$steps])), 
                       na.rm = TRUE)), 
       xlab = "tube number", ylab = "fraction",
       main = paste("pH = ", x$ph),
       type = "h", col = "blue", lwd = 6)
  points(x = 0, y = 1 - sum(x$fractionA[, x$steps]), 
         type = "h", col = "blue", lwd = 6)
  points(x$fractionB[, x$steps], type = "h", col = "green", lwd = 4)
  points(x = 0, y = 1 - sum(x$fractionB[, x$steps]), 
         type = "h", col = "green", lwd = 4)
  points(x$fractionC[, x$steps], type = "h", col = "red", lwd = 2)
  points(x = 0, y = 1 - sum(x$fractionC[, x$steps]),
         type = "h", col = "red", lwd = 2)
  legend(x = "topright", legend = c("A", "B", "C"), 
         lwd = c(6, 4, 2), col = c("blue", "green", "red"), bty = "n")
}