# -*- coding: utf-8 -*-
"""
Created on Tue Nov 21 15:38:23 2023

@author: KAPOTA
"""

# 1- importation des bibliothèques
import numpy as np
import scipy.optimize as op
import matplotlib.pyplot as plt

# 2- Définition des paramètres du problème
DrH0 = (2*-241.8-4*-92.8)*10**3 #J.mol-1
DrS0 = (2*223.1+2*188.7-4*186.9-205.1) #J/K/mol
T0 = 298.00 #K
TF= 2000 #K
R= 8.314 #J/K/mol
P0=1 #bar
P=1 #bar
N=100

# 3- Fonction DrG0 dans l'approx d'Ellingham
def DrG0(T):
    return DrH0 - T*DrS0

# 3- Fonction K dans l'approx d'Ellingham
def K(T) :
    return np.exp(-DrG0(T)/(R*T))

# 3- Fonction à annuler si Q=num/denom=K alors f(T,tau)=num - K*denom
def f(T,P,tau):
    return tau**4*(5 - tau)*P0-K(T)*4**2*(1-tau)**5*P

# 3-définir la fonction dichotomie 
def dichoto(min,max,ff,precision):
    while (max-min) > precision:
        m=(max+min)/2
        if ff(min)*ff(m) < 0 :
            max=m
        else :
            min=m
    return m

# 3-liste des abscisses :
T_abs=np.linspace(T0,TF,N)

# 3- calcul de la liste de tau par la fonction newton
# utilisation de op.newton(g,precision) avec g une fonction à 1 valeur
# on va stocker les valeurs dans la liste tau_newton
tau_newton=[]
for i in range(len(T_abs)):
    def g(tau):
        return f(T_abs[i],P,tau)
    tau_newton.append(op.newton(g,0))
    
# 3- calcul de la liste de tau par la fonction bisect
# utilisation de op.bisect(g,min,max) avec g une fonction à 1 valeur
# on va stocker les valeurs dans la liste tau_bisect
tau_bisect=[]
for i in range(len(T_abs)):
    def h(tau):
        return f(T_abs[i],P,tau)
    tau_bisect.append(op.bisect(h,0,1))

# 3- calcul de la liste de tau par la fonction dichoto
# utilisation de op.bisect(g,min,max) avec g une fonction à 1 valeur
# on va stocker les valeurs dans la liste tau_dichoto
tau_dichoto=[]
for i in range(len(T_abs)):
    def g(tau):
        return f(T_abs[i],P,tau)
    tau_dichoto.append(dichoto(0,1,g,0.0001))
    
# 3- calcul de la liste de tau par la fonction fsolve
# utilisation de op.fsolve(g,min,max) avec g une fonction à 1 valeur
# on va stocker les valeurs dans la liste tau_fsolve
tau_fsolve=[]
for i in range(len(T_abs)):
    def l(tau):
        return f(T_abs[i],P,tau)
    tau_fsolve.append(op.fsolve(l,0.000001))


#4-tracé de la courbe
plt.plot(T_abs,tau_newton,'-r',label='par Newton')
plt.xlabel('T (K)')
plt.ylabel('taux davancement (mol)' )
plt.title('évolution de tau avec T pour l équilibre de Deacon méthode Newton')
plt.grid()
plt.show()

plt.plot(T_abs,tau_bisect,'-b',label='par Bisect')
plt.xlabel('T (K)')
plt.ylabel('taux davancement (mol)' )
plt.title('évolution de tau avec T pour l équilibre de Deacon méthode Bisect')
plt.grid()
plt.show()

plt.plot(T_abs,tau_dichoto,'-g',label='par dichoto')
plt.xlabel('T (K)')
plt.ylabel('taux davancement (mol)' )
plt.title('évolution de tau avec T pour l équilibre de Deacon méthode Dichoto')
plt.grid()
plt.show()

plt.plot(T_abs,tau_fsolve,'-y',label='par Newton')
plt.xlabel('T (K)')
plt.ylabel('taux davancement (mol)' )
plt.title('évolution de tau avec T pour l équilibre de Deacon méthode fsolve')
plt.grid()
plt.show()