#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 29 23:01:37 2017

@author: Alain
"""

import matplotlib.pyplot as plt


""" courbes de répartition diacide """

def pct_H2A(pH, pK1, pK2):
    return 100/(1 + 10**(pH-pK1)+10**(2*pH-pK1-pK2))
    
def pct_HA(pH, pK1, pK2):
    return 100/(10**(pK1-pH) + 1 + 10**(pH-pK2))
    
def pct_A(pH, pK1, pK2):
    return 100/(10**(pK1 + pK2 - 2*pH) + 10**(pK2-pH) + 1 )
 
H2A = []
HA = []
A = []
PH = []
pKa1 = 4
pKa2 = 8
pH = 0
while pH < 14 : 
    PH.append(pH)
    H2A.append(pct_H2A(pH, pKa1, pKa2))
    HA.append(pct_HA(pH, pKa1, pKa2))
    A.append(pct_A(pH, pKa1, pKa2))
    pH = pH + 0.01
    
plt.plot(PH, H2A, 'r', label = "H2A")
plt.plot(PH, HA, 'g', label = "HA-")
plt.plot(PH, A, 'b', label = "A2-")
plt.ylabel("Pourcentage")
plt.xlabel("pH")
plt.show()

""" Courbes de répartition : Fe(s)/Fe2+/Fe3+ """

N0 = 0.01
E1 = - 0.44
E2 = 0.77

def nFe2(E):
    return N0/(1+10**((E-E2)/0.06))
def nFe3(E):
    return N0/(1+10**(-(E-E2)/0.06))
def nFe2_Fe(E):
    return 10**((E-E1)/0.03)
def nFe3_Fe(E):
    return nFe2_Fe(E)*10**((E-E2)/0.06)  
def nFe(E):
    return N0 - (nFe2_Fe(E)+nFe3_Fe(E))

E=[]
FE=[]
FE2=[]
FE3=[]

pot = -0.6
while pot < 1 :
    E.append(pot)
    if nFe(pot)>=0:
        FE.append(nFe(pot))
        FE2.append(nFe2_Fe(pot))
        FE3.append(nFe3_Fe(pot))
    else :
        FE.append(0)
        FE2.append(nFe2(pot))
        FE3.append(nFe3(pot))
    pot = pot + 0.01
    
plt.plot(E, FE2, 'b', label=r"Fe$^{2+}$")
plt.plot(E, FE3, 'r', label=r"Fe$^{3+}$")
plt.plot(E, FE, 'g', label="Fe(s)")
plt.ylabel("Quantité de matière")
plt.xlabel("E(V)")
plt.legend()
plt.show()


"""   courbes de répartition polyacide """

monoacide = ["A", -1, 4]
diacide = ["A", -2, 4, 8]
triacide = ["A", -3, 2, 7, 12]


def c_base(C0, pH, pKa):
    nb_acide = len(pKa)
    denum = 1 
    for i in range(1, nb_acide+1):
        s = 0
        for j in range(len(pKa)-1, len(pKa)-1-i, -1) :
            s = s + pKa[j]
        denum = denum + 10**(s - i*pH)
    c = C0/denum  
    return c
    
def nom(nbH, nomBase, chargeBase):
    name = ""
    if nbH > 0 : 
        name = "H"
        if nbH >= 2 :
            name = name + str(nbH)
    name = name + nomBase
    charge = chargeBase +nbH
    if charge == -1 : 
        name = name + "-"
    elif charge == 1 : 
        name = name + "+"
    elif charge < 0 : 
        name = name + str(abs(charge)) + "-"
    elif charge > 0 : 
        name = name + str(charge) + "+"
    return name
    
    
def plot_repartition(acide):
    pKa = [acide[i] for i in range(2, len(acide))]
    nb_acide = len(acide) - 2
    for i in range (nb_acide + 1):
        PH = []
        PCT = []
        pH = 0
        while pH <= 14 :
            PH.append(pH)
            c = c_base(100, pH, pKa)
            for j in range(i) :
                c = c * 10**(pKa[len(pKa)-j-1]-pH)
            PCT.append(c)
            pH = pH + 0.01
        plt.plot(PH, PCT, label = nom(i, acide[0], acide[1]))
    
plot_repartition(triacide)
plt.ylabel("Pourcentage")
plt.xlabel("pH")
plt.legend()
plt.show()


""" MLL's challenge """

Histidine = ["Histidine", -1, 1.8, 6.1, 9.2]
Glutamique = ["Glutamique", -2, 1.8, 4.5, 9.8]
Leucine = ["Leucine", -1, 2.4, 9.6]


def plot_charge(acide, C0):
    pKa = [acide[i] for i in range(2, len(acide))]
    nb_acide = len(acide) - 2
    PH = []
    Q = []
    pH = 0
    while pH <= 14:
        PH.append(pH)
        q = 0
        for i in range(nb_acide +1):
            c = c_base(C0, pH, pKa)
            for j in range(i) :
                c = c * 10**(pKa[len(pKa)-j-1]-pH)
            q = q + c*(acide[1]+i)
        Q.append(q)
        pH = pH + 0.01
    plt.plot(PH, Q, label = acide[0])
    
plot_charge(Histidine, 1)
plot_charge(Glutamique, 1)
plot_charge(Leucine, 1)
plt.ylabel("Charge")
plt.xlabel("pH")
plt.legend()
plt.show()
    


