from matplotlib import pyplot as plt
from math import *
import time
k10=5
km10=1
k20=2
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")
#P=input("point (absisse) tout les 1/... (sachant que le graph iras de 0 a  nombre de point/le pas") #ancienne methode
J=input("echelle en abscisse , de 0 a  ...")
K=int(deltaT)
X=int(Y)
L=int(J)
P=X/L
D=int(P)
#print(D)


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)))

def suite_euler(n,N):
    A=[0 for i in range(n-1)]
    B=[0 for i in range(n)]
    C=[0 for i in range(n)]
    A.insert(0,1)  # c'est la condition initiales sur le reactif , pour les produits elle est nulle.
    for i in range (1,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
    return(A,B,C)

def suite_kutta2(n,N):
    A=[0 for i in range(n-1)]
    B=[0 for i in range(n)]
    C=[0 for i in range(n)]
    A1=[0 for i in range(n)]
    B1=[0 for i in range(n)]
    C1=[0 for i in range(n)]
    A2=[0 for i in range(n)]
    B2=[0 for i in range(n)]
    C2=[0 for i in range(n)]
    A.insert(0,1)
    for i in range(1,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]
    return(A,B,C)
#A,B,C=suite_euler(X,D)
V=[i/D for i in range(X)]
a=time.clock()
D,E,F=suite_kutta2(X,D)
print(time.clock()-a)
plt.plot(V,D,label="reactif+")
plt.plot(V,E,label="produit1+")
plt.plot(V,F,label="produit2+")
#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()