import numpy as np
import matplotlib.pyplot as plt

#paramètres physiques
rho=8960
C=380
K=400
R=0.02
L=0.3
Tgauche=380
Tdroit=320
T0=293
Text=293

#paramètres de simulation
Nx=30 #discrétisation de la variable d'espace
deltax=L/Nx #pas d'espace
tFin=2000 
Nt=10000 #discrétisation de la variable de temps
deltat=tFin/Nt #pas de temps
tRemarquables=[0,Nt//500,Nt//200,Nt//50,Nt//10,Nt//3,2*Nt//3,4*Nt//5,Nt-1]
# choix d'indices temporels pertinents pour l'illustration graphique


def chaleur1D(h):
    D=K/rho/C
    coeff = D*deltat/deltax**2
    coeff2= 2*h*deltat/R/rho/C
    if  coeff >= 1/2:
        print("Pb de stabilité, réponse sans doute erronée") 
    u=np.zeros((Nx+1,Nt+1))
    u[0]=T0 #initialise toute la 1ère ligne (t=0)
    for j in range(Nt):
        u[0,j+1]=Tgauche # bord imposé
        u[Nx,j+1]=Tdroit # bord imposé
        for i in range(1,Nx):
            der2u=u[i-1,j]-2*u[i,j]+u[i+1,j]
            u[i,j+1]=u[i,j]+coeff*der2u-coeff2*(u[i,j]-Text)
    return u
    
def profils(h,evol):
    X=np.linspace(0,L,Nx+1)
    plt.figure("Quelques profils de températures")
    plt.title("h="+str(h)+", L="+str(L)+", R="+str(R))
    for j in tRemarquables:
        plt.plot(X,evol[:,j],label="t="+str(round(j*deltat)))
    plt.legend()
    
def carteTemp(h,evol):
    tFin2=tFin/10
    Ntbis=Nt//10 #suivi moins long pour la carte de couleurs
    pasdet=Nt//500 #pas plus longs
    plt.figure("Evolution des températures")
    plt.title("h="+str(h)+", L="+str(L)+", R="+str(R))
    valeurs=evol-Text
    selectionValeurs=valeurs[:,:Ntbis:pasdet]
    plt.xticks([0,Ntbis//pasdet-1],[0,tFin2])
    plt.imshow(selectionValeurs)
    
def simulation(h):
    evol=chaleur1D(h)
    profils(h,evol)
    carteTemp(h,evol)
    plt.show()

simulation(400)