import numpy as np
import matplotlib.pyplot as plt
R=8.314

## Fonction cinétique :
tau=1 # minute

def constante_vitesse(T):
    '''renvoie la valeur de la constante de vitesse a la temperature T (en Kelvin). La constante des gaz parfaits R est une variable globale, definie en dehors de la fonction'''
    return np.exp(17-48500/(R*T))  # pour T scalaire ou tableau...

def X_cinetique(T):
    ''' renvoie le taux de conversion pour tau fixe ; tau est variable globale definie hors de la fonction'''
    k=constante_vitesse(T) # calcul de la constante de vitesse
    return tau*k/(1+tau*k) # taux de conversion pour cinetique d ordre 1


## Fonction thermique pour une marche adiabatique :

DrH= - 220000 # enthalpie de reaction exothermique (J/mol)
rho=1000 # g.L-1 masse volumique du fluide entrant
cp=4.0 # J.K-1.g-1 capacite thermique massique du fluide entrant
CAe=2.0 # concentration de A (mol.L-1)
J = - (CAe*DrH)/(rho*cp) # elevation adiabatique : 1/pente de XA=f(T)
Te=290

def X_adia(T):
    ''' renvoie le taux de conversion pour transformation adiabatique. La temperature d'entree Te et l'elevation adiabatique J sont des variables globales'''
    return (T-Te)/J


## Fonction thermique pour une marche avec flux thermique non nul :

N , T0 = 0.3 , 290

def X_thermique(T):
    '''renvoie le taux de conversion pour transformation avec flux thermique . La temperature d'entree Te, l'elevation adiabatique J, la temperature de la paroi exterieure T0 et le rapport N entre tau et tau_c sont des variables globales'''
    return (T-Te)/J + N*(T-T0)/J

## Graphes obtenus pour réacteur adiabatique
Liste_Te=[275,290,310] # liste des temperatures d'entree

plt.close()
T=np.linspace(250,450,100)
for Te in Liste_Te:
    plt.plot(T,X_cinetique(T)-X_adia(T),label=Te)
plt.legend()
plt.title('Différence de X_cin et X_adia pour 3 valeurs de Te')
plt.grid()
plt.show()


## Graphes obtenus pour réacteur non adiabatique
Liste_Te=[275,290,310] # liste des temperatures d'entree

plt.figure()
T=np.linspace(250,450,100)
for Te in Liste_Te:
    plt.plot(T,X_cinetique(T)-X_thermique(T),label=Te)
plt.legend()
plt.title('Différence de X_cin et X_th pour 3 valeurs de Te')
plt.grid()
plt.show()


## Fonction de recherche des racines par dichotomie :

def dicho(f,deb,fin,eps):
    ''' fonction de recherche de racine par dichotomie'''
    g=deb
    d=fin
    assert f(g)*f(d)<=0
    while abs(d-g) > eps:
        m=(g+d)/2 
        if f(m)*f(g)<=0:
            d=m
        else:
            g=m
    return (g+d)/2

## Recherche "automatisée" des racines : 
def cherche_intervalles(diff,deb,fin):
    '''cherche les intervalles de temperatures pour lesquels la fonction diff a une racine : deb designe la temperature initiale et fin designe la temperature finale de l'intervalle d etude'''
    n=50 # nombre de petits intervalles test
    dT=float(fin-deb)/n # pas de temperature (test sur n+1 points)
    LT=[deb] # liste des temperatures avec changement de monotonie
    croissant=diff(deb+dT)-diff(deb)>0 # initialisation du booleen croissant : True si la fonction diff est croissante, False sinon
    for i in range(1,n): # on parcourt tous les petits intervalles dT sauf le premier, deja traite avant la boucle
        T=deb+i*dT # borne gauche de l'intervalle de temperature
        monotonie=diff(T+dT)-diff(T)>0 # monotonie a l abscisse T
        if monotonie!=croissant: # on change de monotonie
            LT.append(T) # ajout de la temperature T
            croissant=monotonie # actualisation de croissant
    LT.append(fin) # on ajoute la derniere temperature
    # verifier la presence d'une racine entre les temperatures de changement de monotonie:
    intervalles=[] # initialisation de la liste des intervalles
    for k in range(len(LT)-1): # parcours de LT par indices
        T=LT[k]
        T_suiv=LT[k+1]
        if diff(T)*diff(T_suiv)<=0: # changement de signe, donc racine presente dans l'intervalle [T,T_suiv]
            intervalles.append((T,T_suiv)) # ajout de l'intervalle : tuple (T,T_suiv)
    return intervalles


def points_fonctionnement(Xa_cin,Xa_therm,Tinf,Tsup):
    '''cherche les points de fonctionnement pour des temperatures entre Tinf et Tsup, et  precise s'il s'agit de points stables ou instables. Xa_cin correspond a la fonction renvoyant le taux de conversion obtenu par le bilan de matiere (cinetique) en fonction de T. Xa_therm correspond a la fonction renvoyant le taux de conversion obtenu par le bilan d energie en fonction de T'''
    def diff(T): # fonction difference pour la dichotomie (depend de Te,J,T0,..)
        return Xa_cin(T)-Xa_therm(T)
    # determination rapide des intervalles de recherche
    intervalles=cherche_intervalles(diff,Tinf,Tsup)
    # recherche precise des solutions
    solutions=[]
    dT=0.2 # pour le test de la stabilite du point (ecarte de dT)
    for tup in intervalles:
        a,b=tup # bornes inferieure et superieure de l'intervalle
        sol=dicho(diff,a,b,1e-2) # temperature de fonctionnement
        x=Xa_therm(sol) # taux de conversion de fonctionnement
        if Xa_cin(sol+dT)>Xa_therm(sol+dT): # point instable 
            stable='point instable'
        else:
            stable='point stable'
        solutions.append(((sol,x),stable)) # ajout du tuple (point de fonctionnement,stabilite du point)
    return solutions

## Application dans le cas de 3 températures d'entrée pour un réacteur adiabatique :
Liste_Te=[275,290,310] # liste des temperatures d'entree

plt.close()
T=np.linspace(250,450,100)
for Te in Liste_Te:
    plt.plot(T,X_cinetique(T)-X_thermique(T),label=Te)
plt.legend()
plt.title('Différence de X_cin et X_th pour 3 valeurs de Te')
plt.grid()
plt.show()

for Te in Liste_Te:
    print('pour Te= ',Te, ' K :')
    resultat=points_fonctionnement(X_cinetique,X_adia,250,450)
    for tup in resultat:
        print(tup[0],' : ',tup[1])
    print() # pour sauter une ligne entre differentes valeurs de Te

#On obtient les résultats suivants (cohérents avec les courbes) :
# pour Te= 275 K : (276.84, 0.0167) : point stable
# pour Te= 290 K : (297.58, 0.0689) : point stable (338.21, 0.4383) : point instable (385.05, 0.8641) : point stable
# pour Te= 310 K : (414.38, 0.9489) : point stable

## Application dans le cas de 3 températures d'entrée pour un réacteur non adiabatique :
Liste_Te=[275,290,310] # liste des temperatures d'entree

for Te in Liste_Te:
    print('pour Te= ',Te, ' K :')
    resultat=points_fonctionnement(X_cinetique,X_thermique,250,450)
    for tup in resultat:
        print(tup[0],' : ',tup[1])
    print() # pour sauter une ligne entre differentes valeurs de Te

#On obtient les résultats suivants (cohérents avec les courbes) :
# pour Te= 275 K : pour Te=  275  K : (280.29, 0.0216)  :  point stable
# pour Te= 290 K : (294.95, 0.0586) : point stable
# pour Te= 310 K : (372.34, 0.79)  :  point stable : point stable