#Cinétique de la formation d'une couche de GaAs

import numpy as np
import matplotlib.pyplot as plt

#1ère partie : Evolution temporelle de [AsH3]= CA
#Données
CA0= 1                 #concentration initiale en AsH3    mol/L
kd1 = 0.01             # constante de vitesse à T = 580K
Ea= 184                # energie d'activation en kJ/mol
R=8.314                #constante des gaz parfaits

#Calcul de kd à T = 750
kd2 = kd1*np.exp((-Ea*1000/R)*(1/750-1/580))
print("valeur de kd à 750K =",kd2,"s-1")

# Détermination de CA(t) à 750 K  et tracé

# 1) Méthode D'Euler
def Fe(C,t):
    return (-kd2*C)
def euler(Fe,t0,tf,CA0,p):
    h=(tf-t0)/(p-1)
    LC = [CA0]
    LT = [0]
    c=CA0
    d=t0
    for k in range(p):
        c=c+h*Fe(c,d)
        d=d+h
        LC.append(c)
        LT.append(d)
    return (LT,LC)
dates,Concentrations = euler(Fe,0,0.1,CA0,10)

# 2) Fonction odeint
from scipy.integrate import odeint
n=200       #nombre de valeurs de la variable t dans l'intervalle retenu
T=np.linspace ( 0,0.1 ,n)
def F(CA,t):
    return (-kd2*CA)
solCA = odeint ( F, [CA0],T)

# 3) réolution analytique
def CA3(t):
    return (CA0*np.exp(-kd2*t))
CAan=CA3(T)

plt.figure(0)
plt.plot(dates,Concentrations,'xk',label='CA(t) Euler')
plt.plot(T,solCA,'ob',label='CA(t)')
plt.plot(T,CAan,'+ r',label='CA(T) analytique')
plt.legend()
plt.xlabel("t(s)")
plt.ylabel ("CA(mol/l)")
plt.show()


# 2ème partie : épaisseur en fonction du temps
#Données
CT0= 0.01                  #concentration initiale en Ga(CH3)3   mol/L
CM0=0                     #concentration initiale en  Ga(CH3)    mol/L
k1 = 5e-2                  # constante de vitesse à T = 750K
k2 =5e-4                   # energie d'activation en kJ/mol
alpha=10**4
mv= 5.3e6                 # masse volumique  en g/m3
M= 145                    #masse molaire de GaAs en g/mol

n=200       #nombre de valeurs de la variable t dans l'intervalle retenu
T2=np.linspace ( 0,500 ,n)

# Système d'équations différentielles
def F2(X,t):
    CT,CM,e =X[0],X[1],X[2]
    return (-k1*CT,k1*CT-k2*CM,k2*1000*CM/alpha/mv)

sol=odeint(F2,[CT0,CM0,0],T2)

CT=sol[:,0]
CM=sol[:,1]
e=sol[:,2]

plt.figure(1)
plt.plot(T2,e,'xr',label='épaisseur')
plt.legend()
plt.xlabel("t(s)")
plt.ylabel ("e(m)")
plt.show()

#taux de croissance
def tau(C):
    return (k2*1000*C/alpha/mv)
CM2=[]
for i in range (len(T2)):
    CM2.append (CM[i])
Y=[]
for i in range (len(T2)):
    Y.append (tau(CM2[i]))


plt.figure(2)
plt.plot(T2,Y,'xg',label='taux de croissance')
plt.legend()
plt.xlabel("t(s)")
plt.ylabel ("de/dt")
plt.show()

# Instant t pour lequel le taux est maximal
#résolution analytique
tmax = (1/(k2-k1))*np.log(k2/k1)
print ("tmax analytique = ",tmax,"s")

# A partir de la résolution numérique
print ("tmax numerique =",T2[np.argmax(Y)])