import numpy as np
import matplotlib.pyplot as plt

## Paramètres physiques

rho = 8960         # Masse volumique (km/m3)
C = 380            # Capacite thermique (J/K/kg)
K = 400            # Conductivite thermique (W/K/m)
H = 100            # Coefficient conducto-convectif (W/K/m2)
R = 0.2            # Rayon du cylindre (m)
L = 0.3            # Longueur du cylindre (m)
Pv = 0.1           # Puissance volumique (W/m3)
T0 = 293           # Temperature initiale (K)
Text = 293         # Temperature exterieure (K)
Tg = 380           # Temperature du thermostat de gauche (K)
Td = 320           # Temperature du thermostat de droite (K)

## Paramètres de la simulation

Nx = 30            # Discrétisation de la variable d'espace
deltax = L/Nx      # Pas d'espace
tFin = 2000        # Temps de fin de simulation (s)
Nt = 10000         # Discrétisation de la variable de temps
deltat = tFin/Nt   # Pas de temps

ttrace = [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

## Resolution de l'equation de la diffusion

def chaleur1D():
    D=K/rho/C
    AA = D*deltat/deltax**2
    BB = 2*H*deltat/R/rho/C
    CC = Pv*deltat
    # verification du critere de convergence

    u=np.zeros((Nt+1,Nx+1))

    ... #initialisation du profil a t=0 (a completer)
    for i in range(Nt):
        ... # bord imposé (a completer)
        ... # bord imposé (a completer)

        ... # boucle principale en j (a completer)
    return u

## Trace des resultats

def profils(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 i in ttrace:
        plt.plot(...) # a completer ; prévoir une légende
    plt.legend()

def carteTemp(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=np.transpose(evol-Text)
    selectionValeurs=valeurs[:,:Ntbis:pasdet]
    plt.xticks([0,Ntbis//pasdet-1],[0,tFin2])
    plt.imshow(selectionValeurs)

## Fonction de simulation

def simulation():
    evol=chaleur1D()
    profils(evol)
    carteTemp(evol)
    plt.show()

... # lancement de la simulation (a completer)