import numpy as np
import matplotlib.pyplot as plt

## Parametres du système et de la réaction A + B = C (A : anhydride ; B : eau ; C : acide acétique)
# Quantités de matière initiales en mol
nA0 = 1
nB0 = 10
# Température initiale en K
T0 = 310
# Enthlapie standard de réaction en J/mol
deltarH0 = -56e3
# Capacité thermique du réacteur en J/K
Créac = 250 
# Capacités thermiques molaires en J/(mol.K)
CmA = 189.7 
CmB = 75.4
CmC = 119.3
# Capacité thermique du système à t = 0 en J/K
Ctot0=Créac+nA0*CmA+nB0*CmB
print(Ctot0)

# Fonction capacité thermique du système en J/K
def Ctot(ksi):
    return(Ctot0+ksi*(2*CmC-CmA-CmB))

# Facteur pré-exponentiel en 1/min et énergie d'activitation en J/mol
KA = 1.68e7
EA = 50.5e3
# Constante des gaz parfaits
R = 8.314

# Loi d'Arrhénius
def k(T):
    return(KA*np.exp(-EA/(R*T)))

## Résolution numérique
# Nombre d’iterations, instant final et pas de temps en min
n = 1001
tf = 1.5/k(T0)    # on identifie un temps caractéristique associé à la cinétique de la réaction via la constante de vitesse à T0
dt = tf/(n-1)

# Listes des instants
t = [i*dt for i in range(n)]
# Création de listes pour l'avancement et pour la température
av = [0 for i in range(n)]
Temp = [0 for i in range(n)]

# Initialisation
av[0] = 0
Temp[0] = T0
# Resolution par la methode d’Euler
for i in range(n-1):
    dav = dt*k(Temp[i])*(nA0-av[i])
    av[i+1] = av[i] + dav
    Temp[i+1] = Temp[i]-deltarH0*dav/Ctot(av[i])
# Capacité thermique du système à t = tf en J/K
print(Ctot(av[-1]))   # à comparer à Ctot(0)

# Representation graphique T(t)
plt.figure(1)
plt.plot(t,Temp)
plt.xlabel("t (min)")
plt.ylabel("T (K)")
plt.grid()
plt.show()

# Representations graphiques avancement(t) et T(t) dans une même fenêtre graphique 
"""plt.subplot(1,2,1)
plt.plot(t,av)        # on vérifie que l'avancement maximal (=nA0) est atteint
plt.xlabel("temps en min")
plt.ylabel("avancement (mol)")
plt.grid()
plt.subplot(1,2,2)
plt.plot(t,Temp)
plt.xlabel("temps en min")
plt.ylabel("température (K)")
plt.grid()
plt.show()"""