From 560c6a6a6285a847aa2be1aacd1784dd268ebeff Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Mon, 6 Jul 2026 11:59:34 +0200 Subject: [PATCH] Aggiunge animazione 3D isometrica del barattolo con colormap di temperatura Estrude il campo x-z della superficie esterna lungo la circonferenza tramite la stessa gaussiana usata per l'attenuazione del flusso, con shading per rendere visibile la curvatura del cilindro anche nelle zone a temperatura ambiente. --- plot_animazione_3d.py | 150 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 150 insertions(+) create mode 100644 plot_animazione_3d.py diff --git a/plot_animazione_3d.py b/plot_animazione_3d.py new file mode 100644 index 0000000..9207912 --- /dev/null +++ b/plot_animazione_3d.py @@ -0,0 +1,150 @@ +# Animazione 3D del barattolo in vista isometrica: colormap della temperatura +# sulla superficie esterna. +# +# Il modello risolve solo la sezione x-z (vedi CLAUDE.md): la coordinata +# circonferenziale y non è simulata, è collassata in un'attenuazione +# gaussiana del flusso. Per estrudere il campo attorno alla circonferenza si +# riusa la stessa gaussiana: la sovratemperatura rispetto al piano y=0 (dove +# si trova il sensore) viene scalata in funzione della distanza +# circonferenziale dal percorso delle sorgenti. È quindi una ricostruzione +# visiva, non un calcolo di diffusione in y. + +import random + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import cm +from matplotlib.animation import FuncAnimation, PillowWriter +from pathlib import Path + +from config import FASCETTA, SIMULAZIONE +from plot_animazione import DT_FRAME_S as DT_FRAME_S_BASE +from plot_animazione import T_FINE_ANIMAZIONE_S, T_INIZIO_ANIMAZIONE_S, simula_campi +from simulate import configurazione_randomizzata + +# Tempo simulato tra un fotogramma e il successivo dell'animazione 3D. Più +# largo di DT_FRAME_S_BASE perché ricreare una superficie 3D a ogni +# fotogramma è più costoso della semplice imshow 2D. +DT_FRAME_S = 0.15 + +# Millisecondi tra i fotogrammi in riproduzione. +INTERVALLO_RIPRODUZIONE_MS = 60 + +# Numero di punti lungo la circonferenza per l'estrusione della superficie. +N_THETA = 72 + + +def attenuazione_circonferenziale( + y_m: np.ndarray, offset_y_m: float, sigma_m: float +) -> np.ndarray: + # Fattore che scala la sovratemperatura (T - T_ambiente) del piano y=0 + # in funzione della distanza circonferenziale y dal piano stesso, + # normalizzato in modo da valere 1 in y=0. + sigma = max(sigma_m, 1e-9) + esponente = -0.5 * ((y_m - offset_y_m) ** 2 - offset_y_m ** 2) / (sigma * sigma) + return np.exp(esponente) + + +def main() -> None: + rng = random.Random(SIMULAZIONE["seed"]) + cfg_run = configurazione_randomizzata(1, rng) + dati = simula_campi(cfg_run) + + tempi = dati["tempi"] + indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S)) + passo = max(1, round(DT_FRAME_S / DT_FRAME_S_BASE)) + indici_frame = list(range(indice_inizio, len(tempi), passo)) + + T_ambiente = cfg_run["aria"]["temperatura_ambiente_C"] + sorgente = cfg_run["sorgente"] + sigma_m = sorgente["sigma_punto_m"] + offset_y_m = sorgente["offset_y_percorso_m"] + numero_sorgenti = sorgente.get("numero_sorgenti", 1) + distanza_m = sorgente.get("distanza_sorgenti_m", 0.0) + + raggio_m = (FASCETTA["diametro_mm"] / 1000.0) / 2.0 + x_centri_m = dati["x_centri_mm"] / 1000.0 + lunghezza_m = dati["lunghezza_mm"] / 1000.0 + + theta = np.linspace(-np.pi, np.pi, N_THETA) + y_circ_m = theta * raggio_m + attenuazione = attenuazione_circonferenziale(y_circ_m, offset_y_m, sigma_m) + + Xm, Thetam = np.meshgrid(x_centri_m, theta) + Ym = raggio_m * np.sin(Thetam) + Zm = raggio_m * np.cos(Thetam) + + T_max = max(c[:, 0].max() for c in dati["campi"]) + # vmin più basso della temperatura ambiente reale: altrimenti le zone + # fredde cadrebbero sul nero puro di "inferno" e, essendo lo shading + # moltiplicativo, nessuna illuminazione basterebbe a renderle visibili. + norm = matplotlib.colors.Normalize( + vmin=T_ambiente - 0.4 * (T_max - T_ambiente), vmax=T_max + ) + cmap = matplotlib.colormaps["inferno"] + lightsource = matplotlib.colors.LightSource(azdeg=315, altdeg=45) + + fig = plt.figure(figsize=(9, 7)) + ax = fig.add_subplot(projection="3d") + ax.view_init(elev=35.264, azim=45) + ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m)) + ax.set_xlabel("x [m]") + ax.set_axis_off() + + # La colorbar mostra il range reale delle temperature: il norm esteso + # verso il basso serve solo a schiarire il colore di base della + # superficie fredda, non deve comparire nella scala mostrata all'utente. + norm_colorbar = matplotlib.colors.Normalize(vmin=T_ambiente, vmax=T_max) + mappabile = cm.ScalarMappable(cmap=cmap, norm=norm_colorbar) + mappabile.set_array([]) + barra = fig.colorbar(mappabile, ax=ax, shrink=0.6, pad=0.05) + barra.set_label("T [°C]") + + def disegna_frame(k: int): + ax.cla() + ax.view_init(elev=35.264, azim=45) + ax.set_box_aspect((lunghezza_m, 2 * raggio_m, 2 * raggio_m)) + ax.set_axis_off() + + T_lato_esterno = dati["campi"][k][:, 0] + T_superficie = T_ambiente + (T_lato_esterno[None, :] - T_ambiente) * attenuazione[:, None] + colori = cmap(norm(T_superficie)) + + ax.plot_surface( + Xm, Ym, Zm, facecolors=colori, rstride=1, cstride=1, + antialiased=False, shade=True, lightsource=lightsource, linewidth=0, + ) + + x_sorgenti_m = dati["x_riferimenti"][k] + np.arange(numero_sorgenti) * distanza_m + visibili = (x_sorgenti_m >= 0.0) & (x_sorgenti_m <= lunghezza_m) + if visibili.any(): + theta_sorgente = offset_y_m / raggio_m + ax.scatter( + x_sorgenti_m[visibili], + np.full(int(visibili.sum()), raggio_m * np.sin(theta_sorgente) * 1.05), + np.full(int(visibili.sum()), raggio_m * np.cos(theta_sorgente) * 1.05), + color="cyan", s=25, depthshade=False, + ) + + ax.set_title(f"t = {tempi[k]:.3f} s") + return () + + animazione = FuncAnimation( + fig, disegna_frame, frames=indici_frame, + interval=INTERVALLO_RIPRODUZIONE_MS, blit=False, + ) + + if matplotlib.get_backend().lower() == "agg": + percorso = Path("dataset") / "animazione_3d.gif" + animazione.save( + percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS) + ) + print(f"Backend non interattivo: animazione salvata in {percorso}") + return + + plt.show() + + +if __name__ == "__main__": + main()