from matplotlib import pyplot as plt


P=200

def dicho(xmin,xmax,x,p,f,a):
    m=0
    while xmax-xmin>p:
        m=(xmin+xmax)/2
        if x==f(m,a):
            return(m)
        elif (f(m,a)-x)*(f(xmin,a) - x) > 0 :
            xmin=m
        else:
            xmax=m

    return(m)


def f(x,a):
    b=P*P*0.60914*10**(-5)-((2*x)**2)*((1+a-2*x)**2)/((a-x)*(1-3*x)**3)
    return (b)

def tau(a):
    if a<(1/3):
        ksi=dicho(0,a,0,0.000001,f,a)
        tau=3*ksi

    else:
        ksi=dicho(0,1/3,0,0.000001,f,a)
        tau=3*ksi
    return(tau)




A=[0.0001+i*0.0001 for i in range(30000)]
B=[tau(a) for a in A]
plt.plot(A,B)
plt.title("fraction de NH3 à l'équilibre en fonction de la composition initiale; T=723K P=200bars")
plt.xlabel("proportion de N2/H2")
plt.ylabel("Tau de conversion H2")
plt.show()