#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov  8 13:24:20 2017

@author: Alain
"""

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 11 17:14:09 2014

@author: Alain
"""
from math import sqrt, log, exp
import matplotlib.pyplot as plt

"""
routines de régression linéaire
"""
   
"""
dichotomie
"""
def dichotomie(f, xmin, xmax, delta):
    while xmax - xmin > delta:
        if f(xmin)*f((xmin+xmax)/2) < 0:
            xmax = (xmin+xmax)/2
        else:
            xmin = (xmin+xmax)/2
    return (xmin+xmax)/2

"""
diagramme binaire
"""
R = 8.314
TA = 384    #toluène
LA = 38000
TB = 373    #eau
LB = 42000

def rosee(x,Teb, Lvap):
    return Teb/(1 - (R*Teb*log(x)/Lvap))

def fonction(x):
    return rosee(x, TB, LB)-rosee(1-x, TA, LA) # attention au 1-x !!!

XH = dichotomie(fonction, 0.1, 0.9, 0.001)
TH = rosee(XH,TB,LB)

print('Hétéroazéotrope : xH =', XH, ' TH = ', TH)  

Xeau = [XH + i *(1 - XH)/100 for i in range(101) ]
Roseau = [rosee(x,TB,LB) for x in Xeau ]
Xtol = [i *XH/100 for i in range(100) ]
Rostol = [rosee(1-x,TA,LA) for x in Xtol ]

plt.plot(Xeau, Roseau, 'r')
plt.plot(Xtol, Rostol, 'b')
plt.plot([0,1], [TH,TH], 'g')
plt.axis([0, 1, 340, 390])
plt.annotate('Gaz', (0.5, 380) )
plt.annotate('Eau(l) + Toluène(l)', (0.4, 345) )
plt.annotate('Eau(l) + Gaz', (0.75, 360) )
plt.annotate('Toluène(l) + Gaz', (0.1, 360) )
plt.show()


"""   SYNTHESE DE L'AMMONIAC """

#variables globales utilisées dans l'ensemble du programme
R = 8.314
n1 = 1.0 # quantité initiale en N2
n2 = 3.0 # quantité initiale en H2
n3 = 0.0 # quantité initiale en NH3
Tmin = 400.0 #température minimale de l'étude
Tmax = 800.0 #température maximale de l'étude
deltaT = 10.0 #intervalle de température entre deux calculs

DrH0 = - 92000 # enthalpie standard de réaction en J/mol
DrS0 = - 200 # entopie standard de réaction en J/K/mol

def A(ksi, T, P):#calcul de l'affinité chimique de la réaction en fonction de l'avancement, de T, de P
	return -(DrH0 - T*DrS0) -R*T*log(pow(n3+2*ksi,2)*pow(n1+n2+n3-2.0*ksi,2.0)/((n1-ksi)*pow(n2-3.0*ksi,3.0)*P*P))

def ksi_eq(mini, maxi, precision, T, P):#calcul de l'avancement à l'équilibre entre mini et maxi avec la precision souhaitée
    while maxi - mini > precision:
        milieu = (maxi + mini)/2
        if A(mini, T, P)*A(milieu, T, P) <0 :
            maxi = milieu
        else:
            mini = milieu
    return (maxi + mini)/2
    
def plot_to(P):#abscisses et ordonnées d'une courbe taux de conversion = f(T) pour P donnée
    X=[]	
    Y=[]
    T= Tmin
    while T <= Tmax :
        to = ksi_eq(0.00001, min(n1, n2/3), 0.001, T, P)
        X.append(T)
        Y.append(to)
        T = T + deltaT
    return X, Y
   
 
for P in [1, 2, 5]:# on trace la courbe pour les différentes valeurs de la pression (en bar)
    X, Y = plot_to(P)# on récupère abscisses et ordonnées
    plt.plot(X, Y, label = "P = "+str(P))# et on les traces
plt.xlabel('Température (K)') ; 
plt.ylabel('Taux de conversion') ; 
plt.legend()
plt.title('Synthèse ammoniac') ; 
plt.show()

""" TITRAGE D'UN ACIDE FAIBLE """

Ke = pow(10,-14)

""" routines utilitaires pour calculer le pH d'un mélange  """
def f(h,C0,V0,Ct,Vt,Ka): # équation à annuler pour déterminer le pH par dichotomie
    return (h - Ke/h)*(V0 + Vt) + Ct*Vt - C0*V0/(1 + h/Ka)
   
def pH(C0,V0,Ct,Vt,Ka): # méthode de dichotomie pour calculer le pH 
    a=0
    b=14
    fa = f(pow(10,-a), C0, V0, Ct, Vt, Ka)
    while b - a > 0.01 :
        c = (a+b)/2
        fc = f(pow(10,-c), C0, V0, Ct, Vt, Ka)
        if (fa *fc < 0):
            b = c
        else:
            a = c
            fa = fc
    return (a+b)/2

""" programme principal """   
V = [] # tableau de valeur pour stocker les volumes
PH = [] # tableau de valeur pour stocker les pH correspondants
v = 0 # volume "courant"
Ca = 0.01 #concentration en acide titré
Va = 100 #volume titré
Cb = 0.1 #concentration en base
pas = 0.01 #pas d'incrémentation du volume
Vmax = 25 # volume maximal de réactif titrant

while v <= Vmax : 
    V.append(v)
    ppH = pH(0.01, 100, 0.1, v, pow(10,-4))
    PH.append(ppH)
    v = v + pas

plt.plot(V, PH)
plt.show()

