import numpy
import math
from matplotlib.pyplot import *
             

def pC0(pka,ph):
    return -math.log10((10**(-ph)-10**(-14+ph))*(1+10**(-ph+pka)))
             

pka=5
ph = numpy.arange(0,7,0.01)
n=len(ph)
pc = numpy.zeros(n)
for k in range(n):
    pc[k] = pC0(pka,ph[k])
plot(pc,ph)
xlabel('pC')
ylabel('pH')
grid()
             

figure(figsize=(12,10))
for pka in range(0,14,2):
    for k in range(n):
        pc[k] = pC0(pka,ph[k])
    plot(pc,ph,label="pKa = %d"%pka)
axis([0,10,0,7])
xlabel('pC')
ylabel('pH')
grid()
legend(loc="lower right")
             

pka=6
pc = 5
def f(x):
    return pC0(pka,x)-pc
import scipy.optimize
ph = scipy.optimize.bisect(f,0,6.99)
             

def base(pka,ph):
    return 1.0/(1+10**(-ph+pka))
ph = numpy.arange(0,14,0.1)
n=len(ph)
b=numpy.zeros(n)
a=numpy.zeros(n)
pka=5
for k in range(n):
    b[k] = 100*base(pka,ph[k])
    a[k]=100-b[k]
figure(figsize=(10,6))
plot(ph,a,label="AH")
plot(ph,b,label="A-")
xlabel("pH")
ylabel("%")
grid()
legend(loc="upper right")
             

def CMetal(pks,pc0,nn,ph):
    if ph<14.0+(pc0-pks)/nn:
        return 10**(-pc0)
    else:
        return 10**(-pks+nn*(14-ph))
                

cM1 = numpy.zeros(n)
cM2 = numpy.zeros(n)
pc0_1=2.0
pc0_2=3.0
pks=20.0
nn=2
for k in range(n):
    cM1[k] = CMetal(pks,pc0_1,nn,ph[k])/10**(-pc0_1)
    cM2[k] = CMetal(pks,pc0_2,nn,ph[k])/10**(-pc0_2)
figure(figsize=(10,6))
plot(ph,cM1,label="pc0=2")
plot(ph,cM2,label="pc0=3")
xlabel('pH')
ylabel('cM/c0')
legend(loc="upper right")
grid()