import numpy as np
import matplotlib.pyplot as plt

""" Paramètres du problème """
L=0.5  # Longueur du système (en m)
D=1e-5  # Coefficient de diffusion thermique de l'acier (en m^2/s)
T1i=273+60  # Température initiale en x=0 (en K)
T1f=273+30  # Température finale en x=0 (en K)
T2=273+20  # Température permanente en x=L (en K)
delta_t=1*3600  # Durée de passage de T1i à T1f (en s)

T=delta_t*1.5  # Durée totale de la résolution (en s)
Nx=100  # Nombre d'intervalles spatiaux pour la discrétisation
Nt=int(T*D*(Nx**2)/(0.35*(L**2)))  # Nombre d'intervalles temporels pour la discrétisation pour non DV Euler
# Nt=1000   # Pour observer la divergence liée à Euler

print(Nt)

""" Résolution numérique """
x_list=np.linspace(0,L,Nx)
t_list=np.linspace(0,T,Nt)
dx=x_list[1]-x_list[0]
dt=t_list[1]-t_list[0]
T_list=np.zeros((Nt,Nx))

# Condition initiale sur la température
for j in range(0,Nx):
    T_list[0,j]=T1i + (T2-T1i)/(Nx-1) * j

#plt.plot(x_list,T_list[0,:]-273)
#plt.title('Condition initiale')
#plt.xlabel('Position (m)')
#plt.ylabel('Température (en °C)')
#plt.grid()
#plt.show()

# Résolution
for i in range(1,Nt):
    
    T_list[i,Nx-1]=T2
    if t_list[i]<delta_t:
        T_list[i,0]=T1i + (T1f-T1i)/delta_t *t_list[i]
    else:
        T_list[i,0]=T1f
    
    for j in range(1,Nx-1):
        T_list[i,j]=T_list[i-1,j]+D*dt/(dx**2)*(T_list[i-1,j+1]-2*T_list[i-1,j]+T_list[i-1,j-1])
    
""" Affichage de la résolution """
i_affichage=np.linspace(0,Nt-1,5,dtype=int)

#print(T_list[10,:])

for i in i_affichage:
    plt.plot(x_list,T_list[i,:]-273,label='t='+str(t_list[i])+' s')
plt.xlabel('Position (en m)')
plt.ylabel('Température (en °C)')
plt.grid()
plt.legend()
plt.show()