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 (mol)
T0 = 298                    # T0 (K)
t0 = 0                      # temps initial (s)
tF = 5000                   # temps final (s)

## grandeurs liées à la thermodynamique
DrH0 = -20000               # DrH° (J/mol)
cpA = 30                    # C°p,A (J/K/mol)
cpB = 20                    # C°p,B (J/K/mol)
C = 250                     # capacité calorifique du réacteur (J/K)

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

## on définit les listes qui vont recueillir les résultats désirés
Temps=[t0]
Temperature=[T0]
Ksi=[0]

## 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))

# 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):
    return -dksi*DrH0/( (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],Ksi[-1],Delta_t)
    dT=d_T(dksi,Ksi[-1])
    Ksi.append(Ksi[-1]+dksi)
    Temperature.append(Temperature[-1]+dT)

# Calcul de la température adiabatique
T_adia=T0-n*DrH0/(n*cpB+C) ; print(T_adia)

# Tracé du graphe
plt.figure(1)
plt.plot(Temps,Temperature,'g',label="T = f(t)")
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 avec $\rm \Delta_rH$° indépendante de T')
plt.legend()
plt.show()