import numpy as np
import matplotlib.pyplot as plt

cB1 = 84.2 # [=] mM
cB2 = 29.4 # [=] mM
pK1 = 6.03
pK2 = 7.57

pH = np.linspace(6,14,100)

t1 = []
t2 = []
t3 = []
for p in pH:
    t1.append(pow(10,-p))
    t2.append(pow(10,pK1-p)*cB1/pow(1+pow(10,pK1-p),2))
    t3.append(pow(10,pK2-p)*cB2/pow(1+pow(10,pK2-p),2))

beta_all = [np.log(10)*(t1[i] + t2[i] + t3[i]) for i in range(len(pH))]
beta_no_t1 = [np.log(10)*(t2[i] + t3[i]) for i in range(len(pH))]
beta_no_t2 = [np.log(10)*(t1[i] + t3[i]) for i in range(len(pH))]
beta_no_t3 = [np.log(10)*(t1[i] + t2[i]) for i in range(len(pH))]
beta_t1_only = [np.log(10)*(t1[i]) for i in range(len(pH))]
beta_t2_only = [np.log(10)*(t2[i]) for i in range(len(pH))]
beta_t3_only = [np.log(10)*(t3[i]) for i in range(len(pH))]

plt.figure()
ax = plt.subplot(231)
ax.plot(pH, beta_all)
ax.set_title('all')
ax = plt.subplot(232)
ax.plot(pH, beta_t1_only)
ax.set_title('t1')
ax = plt.subplot(233)
ax.plot(pH, beta_t2_only)
ax.set_title('t2')
ax = plt.subplot(234)
ax.plot(pH, beta_t3_only)
ax.set_title('t3 only')

ax = plt.subplot(236)
ax.plot(pH, beta_all)
ax.plot(pH, beta_t1_only)
ax.plot(pH, beta_t2_only)
ax.plot(pH, beta_t3_only)
ax.legend(['all','t1','t2','t3 only'])
plt.show()