import matplotlib.pyplot as plt
import numpy as np
V = 0.5
delta_r_H = -150000
E_a = 157000
A = 1.0*10**15
Dv = 3.0/3600
Te = 200+273.15
S = 0.03
h = 80
u = 900
cp = 2.1
R = 8.31
C_B = 6.16
T_ext=293.15

def alpha_thq(T):
    M=-(u*cp/(delta_r_H*C_B))*(T-(Te-h*S*(T-T_ext)/(u*Dv*cp)))
    return M

def alpha_cin(T):
    M=(V/Dv)*A*np.exp(-E_a/(R*T))
    return M/(1+M)


def dicho (f,a,b,eps,K):

    while b-a>eps:
        m=(a+b)/2
        if (f(m)-K)*(f(a)-K)<0 :
            b=m
        else :
            a=m
        m= (a+b)/2
    return m
def f(T):
    return alpha_thq(T)-alpha_cin(T)

x=np.linspace(310,610,200)
A1 = alpha_thq(x)
A2 = alpha_cin(x)
A3 = f(x)
plt.plot(x,A1,"blue",label='alpha_thq')
plt.plot(x,A2,"red",label='alpha_cin')
#plt.plot(x,A3)
plt.legend()
plt.show()
print("T1 =",dicho(alpha_thq,355,375,0.001,0))
print("T2 =",dicho(alpha_thq,550,570,0.001,1))
print("point d'intersection 1 =", dicho(f,355,375,0.001,0))
print("point d'intersection 2 =", dicho(f,455,475,0.001,0))
print("point d'intersection 3=", dicho(f,550,575,0.001,0))


def derivation(f,c,eps):
    return (f(c+eps)-f(c))/eps

K=-(u*cp+h*S/Dv)/(delta_r_H*C_B)
# K=0.0051623376623376625

print("point d'intersection 1 stable =", derivation(alpha_cin,365.5,0.001) < K)
print("point d'intersection 2 stable =",derivation(alpha_cin,461.6,0.001) < K )
print("point d'intersection 3 stable =",derivation(alpha_cin,558,0.001) < K)