#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 21 12:37:31 2022

@author: julie
"""

#Bibliothèques
import numpy as np
import matplotlib.pyplot as plt

#Caractéristiques
L = 0.5 # Longueur de la barre en mètre
lbda = 50 # à compléter, conducticité du fer
h = 10 # à compléter
mu = 7.86*10**3 # à compléter, masse volumique du fer
c = 444 # à compléter, capacité thermique massique du fer
T0 = 20 # à compléter
T1 = 80 # à compléter

N = 50 # à compléter, Nombre de pas de temps pour la discrétisation temporelle
M = 2000# à compléter, Nombre de point pour la discrétisation spatiale
tmax = 5000  # Durée de l'expérience

Delta_t =  tmax/M# pas de temps en seconde (en accord avec la condition de stabilité), écrire l'expression
Delta_x = L/N # Pas spatial en m

print(Delta_t)
# Calcul de la diffusivité :
D = lbda/mu/c#à compléter

alpha = D*Delta_t/Delta_x**2
print('alpha =', alpha, 'M =', M)

# Tabelau des dates et positions
dates = np.linspace(0,tmax,M+1) # Le tableau des dates
absc = np.linspace(0,L,N+1) #tableau des abscisses

#Matrice des températures 
T = np.zeros((M+1,N+1)) #tableau des écarts de température initiaux

# test qui est qui, pour N et M petits et différents, print('T = ',T)


'Condtion aux limites : bord gauche Tg = T0, bord droit Td=T1'
T[:,0] = T0
T[:,N] = T1

'Condition initiale '
T[0,:] = T0

# test CI et CL pour N et M petits : print(T)

'Boucle temporelle'
for k in range(M) :
    for i in range(1,N) :
        T[k+1,i] = T[k,i] + alpha*(T[k,i+1] + T[k,i-1] - 2*T[k,i])

#tracé toutes les 500 s
    if (k==1 or k==100 or (k*Delta_t)%500 == 0 ) :
        plt.plot(absc,T[k,:], label='t(s)= '+ str(k * Delta_t))

plt.plot(absc,T[M,:], label='t(s)= '+ str(M * Delta_t))

plt.xlabel('x (m)')
plt.ylabel('T (°C)',rotation=0)
plt.title('Diffusion thermique ')
plt.legend()
plt.grid()
plt.show()