################################
#PARAMETRES POUR UTILISATION NORMALE
###################################
fen_abs=6.5
fen_ord=4.5
taille_cara=10
larg_trait=1
taille_mark=4

########################
#PARAMETRES POUR LIVRE DUNOD
##############################
fen_abs=18  #18 
fen_ord=12  #12
taille_cara=30   #30
larg_trait=4      #4
taille_mark=15    #15



import numpy as np               #bibliothèque numpy renommée np
import matplotlib.pyplot as plt  #module matplotlib.pyplot renommé plt


ItMax=5000    #nombre maximal d'itérations
epaiss=40e-2  #épaisseur du mur
conduc=1.65   #conductivité thermique du mur
rho=2150      #masse volumique du mur
cp=1000       #capacité thermique massique du mur
Tint=20       #température intérieure en °C
Text1=10      #température extérieure 1 en °C
Text2=-10     #température extérieure 2 en °C


N=60   #nombre de points de calcul à l'intérieur du mur (60 points)
Dt=25  #intervalle de temps delta t = 25 s
Dx=epaiss/(N+1) #N+2 points, donc N+1 intervalles
alpha=rho*cp/conduc
r=Dt/(alpha*(Dx**2))   #il faut que r < 1/2
a=(Text1-Tint)/epaiss
b=Tint

x=np.zeros(N+2)  #N+2 valeurs en abscisses
for i in range(N+2):
    x[i]=i*Dx


def calc_norme(V): #calcul de la norme du vecteur V
    somme=0
    for i in range(len(V)):
        somme+=V[i]**2
    #on peut remplacer la boucle par : somme=np.sum(V*V)    
    return np.sqrt(somme)
    
def schema_explicite(ItMax,Tint,Text2,a,b,r,N):
    T_tous_k=np.zeros((N+2, ItMax+1))
    T_tous_k[:, 0]=a*x+b #profil de température pour t = 0
    if r>=1/2:
        print("La condition r<1/2 n'est pas respectée.")  
    nbIter=0 #initialisation du nombre d'itérations

    while nbIter<(ItMax-1):
        k=nbIter+1
        T_tous_k[0, k]=Tint #colonne 0 à Tint       
        T_tous_k[N+1, k]=Text2 #colonne N+1 à Text2
        
        #utilisation de la boucle for : PLUS LENT
#        for i in range(1,N+1): #i varie entre 0 et N
#            T_tous_k[i,k]=T_tous_k[i,k-1]*(1-2*r)+(r*T_tous_k[i+1,k-1])+(r*T_tous_k[i-1,k-1])
        
        #utilisation du slicing : BEAUCOUP PLUS RAPIDE
        A=np.zeros(N)
        A=  -2*T_tous_k[1:N+1,k-1]+  T_tous_k[0:N,k-1] +T_tous_k[2:N+2,k-1]
        T_tous_k[1:N+1,k]=T_tous_k[1:N+1,k-1]+r*A

        nbIter+=1 #on incrémente de 1 la variable nbIter
        if calc_norme(T_tous_k[:, k]-T_tous_k[:, k-1])<1e-3: #test si écart < 10^-3
            return nbIter, T_tous_k #quitte la fonction schema_explicite
                
    return nbIter,T_tous_k #retourne le nombre d’itérations et la matrice T_tous_k




#####################################################
##### PROGRAMME PRINCIPAL SCHEMA EXPLICITE
#####################################################
nbIter, T_tous_k=schema_explicite(ItMax, Tint, Text2, a, b, r, N)
t_rp=int(nbIter*Dt/3600)
print('Régime permanent établi au bout de',t_rp,'heures.')

plt.figure()  #nouvelle fenêtre graphique
for k in range(0, nbIter):
    if k%200==0:    
        plt.plot(x, T_tous_k[:, k])
plt.xlabel('x (en m)')
plt.ylabel('T (en °C)')
plt.title("Schéma explicite avec N = 60")
plt.show()   #affiche la figure à l'écran






print("nombre d'itérations =",nbIter)
temps_car=int(alpha*(epaiss**2)/3600)
print('temps caractéristique =',temps_car,'heures.')

#plt.figure()
#plt.plot(x,T_tous_k_imp[:,0])
#plt.plot(x,T_tous_k_imp[:,200])
#plt.plot(x,T_tous_k_imp[:,400])
#plt.plot(x,T_tous_k_imp[:,600])
#plt.plot(x,T_tous_k_imp[:,800])
#plt.plot(x,T_tous_k_imp[:,1000])
#plt.plot(x,T_tous_k[:,2800])
#for k in range(1000,nbIter):
#    if k%200==0:    
#        plt.plot(x,T_tous_k_imp[:,k])
#plt.show()