# -*- coding: utf-8 -*-
"""
Created on Thu Nov  2 14:43:40 2023

@author: cathe
"""

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

# 2- Définition des paramètres du problème
Ea = 50.5*10**3 #J.mol-1
A = 1.68*10**7 #min-1
n0 = 1.00 #mol
n1=50.00 #mol
Meau = 18 #g/mol
V= n1*Meau*10**-3 #L
C0=n0/V
ksi0 = 0.00 #mol/L
T0 = 298.00 #K
R= 8.314 #J/K/mol
t0 = 0#min temps initial
tf2 = 200#min temps final
N = 200#nombre de pas
Dt = (tf2 - t0)/N #min : pas du temps

ccalo=85.3 #J/K
ceau_mass = 4.18 #J/K/g
DrH=-56*10**3 #J/mol

# 3- fonction k
def k(T) :
    return A*np.exp(-Ea/(R*T)) # attention bien mettre les () et les *
#3- dT
def dT(d_ksi):
    return -DrH*d_ksi*V/(ccalo+n1*Meau*ceau_mass)
#3- T adiabatique
Tadiab=T0-DrH*n0/(ccalo+n1*Meau*ceau_mass)

# 3- d_ksi
def d_ksi(T,ksi,Dt):
    return k(T)*(C0-ksi)*Dt

#3-euler : avec une boucle 
t_euler=[t0]
T_euler=[T0]
ksi_euler=[ksi0]
for i in range(N):
    t_euler.append(t_euler[i]+Dt)
    T_euler.append(T_euler[i]+dT(d_ksi(T_euler[i],ksi_euler[i],Dt)))
    ksi_euler.append(ksi_euler[i]+d_ksi(T_euler[i],ksi_euler[i],Dt))

# 3- intégration en cas isotherme
t_abs=np.linspace(t0,tf2,N+1)
def ksi_int_isotherme(t,T):
    return C0*(1-np.exp(-k(T)*t))


# 4- tracé des courbes
plt.plot(t_euler,T_euler,'-g',label='T adiabatique')
plt.axhline(y=Tadiab,color='r',linestyle='--',label='T_max')
plt.grid() #parce que ça fait joli
plt.title("évolution de T pour une transfo adiabatique")
plt.legend()
plt.xlabel('t (min)')
plt.ylabel('T (K)')
plt.show()

plt.plot(t_euler,ksi_euler,'-b',label='cas adiabatique')
plt.axhline(y=C0,color='g',linestyle=':',label='ksi_max')
plt.plot(t_abs,ksi_int_isotherme(t_abs,T0),'--r',label='cas isotherme')
plt.grid() #parce que ça fait joli
plt.title("évolution de ksi pour une transfo adiabatique")
plt.legend()
plt.xlabel('t (min)')
plt.ylabel('avancement volumique (mol/L)')
plt.show()