import numpy as np
import matplotlib.pyplot as plt
from math import *

## grandeurs liées au problème étudié
n = 1                       # nombre initial de moles de A
T0 = 298                    # T0 (K)
t0 = 0                      # temps initial (s)
tF = 4000                   # temps final (s)

## grandeurs liées à la thermodynamique
DrH0 = -20000               # DrH° à la température T0 (J/mol)
cpA = 30                    # Cp°A (J/K/mol)
cpB = 20                    # Cp°B (J/K/mol)
C = 250                     # capacité calorifique de l'enceinte (J/K)

## grandeurs liées à la cinétique
R = 8.314                       # constante des gaz parfaits en uSI
Ea = 20000                      # énergie d'activation en J/mol
k0 = 6e-4                       # constante cinétique à la température T0 (en s-1)

## on définit les listes qui vont recueillir les résultats désirés
Temps=[t0]
Temperature_1=[T0]          # Température dans le cadre de l'approximation d'Ellingham
Ksi_1=[0]                   # avancement dans le cadre de l'approximation d'Ellingham
Temperature_2=[T0]          # Température sans faire l'approximation d'Ellingham
Ksi_2=[0]                   # avancement sans faire l'approximation d'Ellingham

## on définit les paramètres utiles pour la résolution numérique
Npoints = 10000              # nombre d'itérations pour couvrir l'intervalle de temps [t0,tF]

# valeur de la constante de vitesse en fonction de la température
def k(T):
    return k0*np.exp(-Ea/R*(1/T-1/T0))

# valeur de l'enthalpie standard de réaction en fonction de la température
def DrH0_T(T):
    return DrH0 + (cpB-cpA)*(T-T0)

# variation de ksi pendant l'intervalle de temps dt
def d_ksi(T,ksi,dt):
    return k(T)*(n-ksi)*dt

# variation de la température T pendant l'intervalle de temps dt
def d_T(dksi,ksi,T):
    return -dksi*DrH0_T(T)/( (n-ksi-dksi)*cpA + (ksi+dksi)*cpB + C)

# intervalle de temps entre deux estimations de la température
Delta_t = (tF-t0)/Npoints

# Génération des listes t(temps), ksi et T(température)
for i in range(Npoints):
    Temps.append(Temps[-1]+Delta_t)
    dksi=d_ksi(Temperature_1[-1],Ksi_1[-1],Delta_t)
    dT=d_T(dksi,Ksi_1[-1],T0)
    Ksi_1.append(Ksi_1[-1]+dksi)
    Temperature_1.append(Temperature_1[-1]+dT)
    dksi=d_ksi(Temperature_2[-1],Ksi_2[-1],Delta_t)
    dT=d_T(dksi,Ksi_2[-1],Temperature_2[-1])
    Ksi_2.append(Ksi_2[-1]+dksi)
    Temperature_2.append(Temperature_2[-1]+dT)

# Calcul de la température adiabatique
T_adia=T0-n*DrH0/(n*cpB+C) ; print("valeur de la température adiabatique (en k) :",T_adia)

# Tracé du graphe
plt.figure(1)
plt.plot(Temps,Temperature_1,'g',label="T = f(t) avec l'approximation d'Ellingham")
plt.plot(Temps,Temperature_2,'b',label="T = f(t) sans l'approximation d'Ellingham")
plt.plot([t0,tF],[T_adia,T_adia],'r-',label="Température de réaction adiabatique")
plt.xlabel(r'Temps t (s)',fontsize=14)
plt.ylabel(r'Température T (K)',fontsize=14)
plt.title(r'Évolution temporelle de la température en réacteur adiabatique')
plt.legend()
plt.show()