simulatore-induttori/plot_animazione_3d.py

# 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()