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

### Données numériques

# Caractéristiques de la réaction
DrH0 = -150e3 # enthalpie standard de réaction en J/mol
Ea = 157e3 # énergie d'activation en J/mol
A = 1e15 # facteur préexponentiel en s-1
R = 8.314 # constante des GP en J/K/mol

# Caractéristiques du réacteur et du fluide en entrée
V = 0.5 # volume en L
Dv = 3/3600 # débit volumique en L/s
tau = V/Dv # temps de passage en s
ro = 1000 # masse volumique de l'eau en g/L
cp = 1.49 # capacité thermique massique de l'eau en J/g/K
Ce = 0.9 # cocentration en entrée du réactif en mol/L
Te = 423 # température en entrée dans le cas d'un réacteur calorifugé, en K
Te_bis = 623 # température en entrée dans le cas d'un réacteur refroidi, en K
T0 = 298 # température du fluide caloriporteur en K
S = 300e-4 # surface d'échange avec le fluide caloriporteur en m2
h = 50 # coefficient conducto-convectif du fluide caloriporteur en W/m2/K

### Expressions du taux de conversion en fonction de la température

# Expression du taux de conversion résultant de l'étude cinétique :
def k(T) :
    return A*np.exp(-Ea/R/T)
def alpha_Cin (T) :
    return tau*k(T)/(1 + tau*k(T))
    
# Expression du taux de conversion résultant du bilan énergétique :
def alpha_Th (T) :
    return - ro*cp*(T - Te)/(Ce*DrH0)
    
### Détermination des points de fonctionnement

# Tracé des 2 courbes afin de visualiser les points de fonctionnement :
T = np.linspace (200, 700, 10000) # Températures
plt.figure(0)
plt.plot (T, alpha_Cin (T), '-', color = 'blue', label ='Cin')
plt.plot (T, alpha_Th (T), '-', color = 'red', label ='Th')
plt.xlabel ('Température T (K)')
plt.ylabel (r'$\alpha$')
plt.title ('Point(s) de fonctionnement')
plt.grid ()
plt.legend ()
plt.xlim (300, 700)
plt.ylim (0, 1.1)

# Détermination des coordonnées des points de fonctionnement et mise en évidence sur le graphe :
def fonctionnement (T) :
    return alpha_Cin(T) - alpha_Th(T)

Ts1 = bisect (fonctionnement,400,440)
as1 = alpha_Cin(Ts1)
Ts2 = bisect (fonctionnement,440,480)
as2 = alpha_Cin(Ts2)
Ts3 = bisect (fonctionnement,500,540)
as3 = alpha_Cin(Ts3)

plt.figure(1)
plt.scatter (Ts1, as1, s = 50, facecolor = 'none', edgecolor = 'black')
plt.scatter (Ts2, as2, s = 50, facecolor = 'none', edgecolor = 'black')
plt.scatter (Ts3, as3, s = 50, facecolor = 'none', edgecolor = 'black')

plt.plot (T, alpha_Cin (T), '-', color = 'blue', label ='Cin')
plt.plot (T, alpha_Th (T), '-', color = 'red', label ='Th')
plt.xlabel ('Température T (K)')
plt.ylabel (r'$\alpha$')
plt.title ('Point(s) de fonctionnement')
plt.grid ()
plt.legend ()
plt.xlim (300, 700)
plt.ylim (0, 1.1)

plt.text (600, 0.9, r'$\alpha^{s3} = $' + '{:.4}'.format (as3))
plt.text (600, 0.82, r'$T^{s3} = $' + '{:.4}'.format (Ts3) + 'K')
plt.text (600, 0.52, r'$\alpha^{s2} = $' + '{:.4}'.format (as2))
plt.text (600, 0.44, r'$T^{s2} = $' + '{:.4}'.format (Ts2) + 'K')
plt.text (600, 0.14, r'$\alpha^{s1} = $' + '{:.4}'.format (as1))
plt.text (600, 0.06, r'$T^{s1} = $' + '{:.4}'.format (Ts1) + 'K')

### Détermination de la stabilité des points de fonctionnement

# Fonction qui renvoie la valeur de la dérivée d'une fonction
def derivee (f, x) :
    h = 0.001
    return (f(x+h) - f(x))/h

# Fonction pour déterminer la stabilité d'un point de fonctionnement : celle-ci doit retourner 'True' pour un PD stable et 'False' pour un point instable
def stabilite (Ts) :
    if derivee (alpha_Cin,Ts) < derivee (alpha_Th,Ts) :
        return True
    else:
        return False

print ('PF1 stable :' + str(stabilite(Ts1)))
print ('PF2 stable :' + str(stabilite(Ts2)))
print ('PF3 stable :' + str(stabilite(Ts3)))

### Modification des caractéristiques du réacteur : refroidissement par un fluide caloriporteur

# Nouvelle expression du taux de conversion résultant du bilan énergétique :
def alpha_Th_bis (T) :
  return (- h*S/Dv*(T - T0) - ro*cp*(T - Te_bis))/Ce/DrH0

# Visualisation du nouveau point de fonctionnement
plt.figure(2)
plt.plot (T, alpha_Cin (T), '-', color = 'blue', label ='Cin')
plt.plot (T, alpha_Th_bis (T), '-', color = 'green', label ='Th_bis')
plt.xlabel ('Température T (K)')
plt.ylabel (r'$\alpha$')
plt.title ('Nouveau point de fonctionnement')
plt.grid ()
plt.legend ()
plt.xlim (300, 700)
plt.ylim (0, 1.1)

# Détermination des coordonnées du nouveau point de fonctionnement :
def fonctionnement_bis (T) :
    return alpha_Cin(T) - alpha_Th_bis(T)

Ts_bis = bisect (fonctionnement_bis,450,550)
as_bis = alpha_Cin(Ts_bis)

print ('Point de fonctionnement pour un RCPA refroidi :')
print ('Température de sortie : ',Ts_bis,' K')
print ('taux de conversion : ',as_bis)

# Fonction pour déterminer la stabilité du nouveau point de fonctionnement
def stabilite_bis (Ts) :
    if derivee (alpha_Cin,Ts) < derivee (alpha_Th_bis,Ts) :
        return True
    else:
        return False

print ('nouveau PF stable :' + str(stabilite_bis(Ts_bis)))