from matplotlib import pyplot as plt
from math import *
import time


k10=5
km10=1
k20=1.7
km20=0.2
R=8.314
Ea1=70000
Eam1=75000
Ea2=80000
Eam2=77000
T0=273


deltaT=input("ecart de temperature par rapport a 293K(ref)")
Y=input("nombre de point")
J=input("echelle en abscisse , de 0 a  ...")
K=int(deltaT)
X=int(Y)
L=float(J)
P=X/L
D=int(P)


k1=k10*exp(-((Ea1)/(R))*((1/(T0+K))-1/(T0)))
km1=km10*exp(-((Eam1)/(R))*((1/(T0+K))-1/(T0)))
k2=k20*exp(-((Ea2)/(R))*((1/(T0+K))-1/(T0)))
km2=km20*exp(-((Eam2)/(R))*((1/(T0+K))-1/(T0)))


A=[0 for i in range(X-1)]
B=[0 for i in range(X)]
C=[0 for i in range(X)]
A1=[0 for i in range(X)]
B1=[0 for i in range(X)]
C1=[0 for i in range(X)]
A2=[0 for i in range(X)]
B2=[0 for i in range(X)]
C2=[0 for i in range(X)]
A.insert(0,1)


def suite_euler(a,n,N):

    for i in range (a,n):
        A[i]=A[i-1]*(1-(k1+k2)/N)+(km1*B[i-1]+km2*C[i-1])/N
        B[i]=B[i-1]*(1-km1/N)+(k1*A[i])/N
        C[i]=C[i-1]*(1-km2/N)+(k2*A[i])/N
    E=abs(1-A[n-1]-B[n-1]-C[n-1])
    S=abs(C[n-1]-C[n-2])
    return(E,S)


def suite_kutta2(a,n,N):

    for i in range(a,n):
        A1[i]=A[i-1]*(1-(k1+k2)/(2*N))+(km1*B[i-1]+km2*C[i-1])/(2*N)
        B1[i]=B[i-1]*(1-km1/(2*N))+(k1*A[i])/(2*N)
        C1[i]=C[i-1]*(1-km2/(2*N))+(k2*A[i])/(2*N)
        A2[i]=A1[i-1]*(-(k1+k2))+(km1*B1[i-1]+km2*C1[i-1])
        B2[i]=B1[i-1]*(-km1)+(k1*A1[i])
        C2[i]=C1[i-1]*(-km2)+(k2*A1[i])
        A[i]=A[i-1]+(1/N)*A2[i]
        B[i]=B[i-1]+(1/N)*B2[i]
        C[i]=C[i-1]+(1/N)*C2[i]
    E=abs(1-A[n-1]-B[n-1]-C[n-1])
    S=abs(C[n-1]-C[n-2])        # on peut aussi regarder avec B ou A
    return(E,S)


def remplissage(n,N):
    E=0
    i=1
    S=1
    while E<0.02 and i<-20+n and S>0.0005:
        E,S=suite_euler(i,i+20,N)
        i=i+20
    if E>0.02:
        i=1
        Q=1
        while  i<-20+n and Q>0.00008:
            M,Q=suite_kutta2(i,i+20,N)
            i=i+20

    suite_euler(i,n,N)
    V=abs(1-A[n-1]-B[n-1]-C[n-1])
    if V>0.02:
        print("pas assez de points")
    return()


a=time.clock()
remplissage(X,D)
print(time.clock()-a)
V=[i/D for i in range(X)]
plt.plot(V,A,label="reactif")
plt.plot(V,B,label="produit1")
plt.plot(V,C,label="produit2")
plt.legend()
plt.axis([0,L,0,1])
plt.show()