# Synthèse de l'eau  - DLB25-4

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import bisect

# Données
R= 8.314

# Grandeurs standard de la réaction
def DG(T):
    return (- 495000+7.8*T*np.log(T) + 33*T+0.01*T**2)
def DH(T):
    return (- 495000 - 7.8*np.log(T) - 0.01*T**2)
def DS(T):
    return (-40.8-0.02*T-7.8*np.log(T) )
def K(T):
    return (np.exp(-DG(T)/R/T))

print ("valeur de la constante d'équilibre à 298K :", K(298))
print ("valeur de la constante d'équilibre à 1500K:",K(1500))
print ("valeur de la constante d'équilibre à 3000K:",K(3000))

# Courbes illustrant la variation des grandeurs standard en fonction de la temperature
XT=np.linspace(1500,3000,100)
YG=DG(XT)
YH=DH(XT)
YK=K(XT)


plt.figure(0)
plt.plot(XT,YK,'+b',label="constante d'équilibre")
plt.xlabel("T(K)")
plt.ylabel("DrH°")
plt.grid()
plt.legend()
plt.show()

#Composition du système à l'équilibre
def Qr(x,n1,n2,n3):
    return ((1-x)*(n3+2*x)**2/(n1-2*x)**2/(n2-x))
def equation (x,n1,n2,n3,T):
    return Qr(x,n1,n2,n3)-K(T)

xe3=bisect(equation,-0.45,0,args=(0.067,0.033,0.9,3000))
print("l'avancement à l'équilibre dans le 3ème  cas  est:",xe3)

#résolutions vaines :
#xe1=bisect(equation,0,0.167,args=(0.333,0.333,0.333,1500))
#xe2=bisect(equation,0,0.025,args=(0.05,0.05,0.9,1500))
#print("l'avancement à l'équilibre dans le 1er cas  est:",xe1)
#print("l'avancement à l'équilibre dans le 2ème  cas  est:",xe2)


#Tracé des fonctions pour visualiser l'intervalle ou se trouve la solution
#1er cas
X1=np.linspace(0.166496,0.1664995,100)
Y1= equation (X1,0.333,0.333,0.333,1500 )
plt.figure(1)
plt.plot(X1,Y1,'+k',label="1er cas ")
plt.plot(X1,0*X1,'--r')
plt.xlabel("x")
plt.ylabel("Qr-K°")
plt.grid()
plt.legend()
plt.show()

#2ème cas
X2=np.linspace(0.02485,0.024998,100)
Y2= equation (X2,0.05,0.05,0.9,1500 )
plt.figure(2)
plt.plot(X2,Y2,'+b',label="2èmecas ")
plt.plot(X2,0*X2,'--r' )
plt.xlabel("x")
plt.ylabel("Qr-K°")
plt.grid()
plt.legend()
plt.show()

#résolution d el'équation dans le sdeux premiers cas
xe1=bisect(equation,0.166496,0.1664995,args=(0.333,0.333,0.333,1500))
xe2=bisect(equation,0.166496,0.1664995,args=(0.05,0.05,0.9,1500))



print("l'avancement à l'équilibre dans le 1er cas  est:",xe1)
print("l'avancement à l'équilibre dans le 2ème  cas  est:",xe2)