#Réaction d'isomérisation en phase gazeuse - Capacité numérique -Thermo1

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import bisect


#A- Etude isotherme en réacteur fermé
#Données
R=8.314
Ea1=46400
Ea2=118400
A1=2.95e7
A2=1.57e18

# Constantes de vitesse en fonction de la température
def k1(T):
    return (A1*np.exp(-Ea1/R/T))
def k2(T):
    return (A2*np.exp(-Ea2/R/T))

## Q2. Composition du système en fonction du temps- T fixée à 298 K
### Taux de transformation en fonction du temps , T = 298 K / résolution analytique
def k(T):
    return (k1(T)+k2(T))

def alpha1(t):
    return (k1(298)*(1-np.exp(-k(298)*t))/k(298))

temps=np.linspace(0,60,500)    #intervalle de temps
Y1=alpha1 (temps)

plt.figure(0)
plt.plot(temps,Y1,'g',label="taux de conversion , T = 298 K ")
plt.grid()
plt.legend()
plt.xlabel("t(s)")
plt.ylabel("alpha")
plt.show()

print (k1(298))
print (k(298))
print (k1(298)/k(298))


## Q3.Influence de la température
#Q3a.
def alpha298(t) :
    return (k1(298)*(1-np.exp(-k(298)*t))/k(298))
def alpha350(t) :
    return (k1(350)*(1-np.exp(-k(350)*t))/k(350))
def alpha400(t) :
    return (k1(400)*(1-np.exp(-k(400)*t))/k(400))

Y298=alpha298 (temps)
Y350=alpha350 (temps)
Y400=alpha400 (temps)


plt.figure(1)
plt.plot(temps,Y298,'g',label="T=298K")
plt.plot(temps,Y350,'b',label="T=350K")
plt.plot(temps,Y400,'r',label="T=400K")
plt.title("influence de la temprature")
plt.grid()
plt.legend()
plt.xlabel("t(s)")
plt.ylabel("alpha")
plt.show()

## Q3d.
## Taux de conversion à l'équilibre en fonction de la température
def alphaeq (T):
    return (k1(T)/k(T))

Xtemp =np.linspace (250,420,500)   # intervalle de températures
Yeq=alphaeq (Xtemp)

plt.figure(2)
plt.plot(Xtemp,Yeq,'g',label=" alpha(T)")
plt.title("taux de conversion en fonction de la température")
plt.grid()
plt.legend()
plt.xlabel("T(K))")
plt.ylabel("alpha equilibre")
plt.show()

## taux de transformation en fonction de T pour une durée fixée
def alpha(T,duree):
    return (k1(T)*(1-np.exp(-k(T)*duree))/k(T))

Y01= alpha(Xtemp,0.1)
Y1=alpha(Xtemp,1)
Y50= alpha(Xtemp,30)

plt.figure(3)
plt.plot(Xtemp,Yeq,'g',label=" alpha equilibre(T)")
plt.plot(Xtemp,Y01,'r',label=" alpha 0,1s")
plt.plot(Xtemp,Y1,'b',label=" alpha  1s")
plt.plot(Xtemp,Y50,'xk',markersize=1,label="alpha  30s")
plt.title("taux de conversion en fonction de la température")
plt.grid()
plt.legend()
plt.xlabel("T(K))")
plt.ylabel("alpha")
plt.show()

##B- Etude de la marche M adiabatique

from scipy.integrate import odeint

#Données
CR= 250
CpA=30
CpB=20



# Résolution du système d'équations différentielles couplées
#expression de la capacité thermique globale ( dénominateur)
def Cth(x):
    return ((1-x)*CpA+x*CpB+CR)

temps2 = np.linspace (0,2,50)

def F(G,t):
    ALPHA,T = G[0],G[1]
    return (k1(T)-k(T)*ALPHA,72000*(k1(T)-k(T)*ALPHA)/Cth(ALPHA))
sol=odeint (F,[0,298],temps2)

ALPHA = sol[:,0]
TEMP =sol[:,1]

plt.figure(4)
plt.plot(temps2,TEMP,'b',label=" T(temps)")
plt.title(" température en fonction du temps ")
plt.grid()
plt.legend()
plt.xlabel("t(s))")
plt.ylabel("T(K)")
plt.show()

# Prise en compte des fuites thermiques
#Données
U = 500
S= 30
Text= 293

def F2(G2,t):
    ALPHA2,T2 = G2[0],G2[1]
    return (k1(T2)-k(T2)*ALPHA2,(72000*(k1(T2)-k(T2)*ALPHA2)+U*S*(Text-T2))/(Cth(ALPHA2)))
sol2=odeint (F2,[0,298],temps2)

ALPHA2 = sol[:,0]
TEMP2 =sol[:,1]

plt.figure(5)
plt.plot(temps2,TEMP2,'b',label=" T(temps)avec fuite")
plt.plot(temps2,TEMP,'--r',label=" T(temps)sans fuite")
plt.title(" température en fonction du temps ")
plt.grid()
plt.legend()
plt.xlabel("t(s))")
plt.ylabel("T(K)")
plt.show()