94 lines
3.0 KiB
Python
94 lines
3.0 KiB
Python
|
import pandas as pd
|
|||
|
|
import numpy as np
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
from scipy.optimize import curve_fit
|
|||
|
|
|
|||
|
|
# --- Dati ---
|
|||
|
|
df = pd.read_csv("data.csv")
|
|||
|
|
df["time_s"] = df["time since start [ms]"] / 1000.0
|
|||
|
|
|
|||
|
|
T_INF = 22.99 # temperatura ambiente [°C]
|
|||
|
|
T_START = 115.0 # inizio finestra di fit [s]
|
|||
|
|
T1 = 115.0 # riferimento 1° esponenziale (fisso)
|
|||
|
|
T2 = 117.5 # riferimento 2° esponenziale (fisso)
|
|||
|
|
W_ZERO_START = 115.9
|
|||
|
|
W_ZERO_END = 117.2
|
|||
|
|
|
|||
|
|
mask = df["time_s"] >= T_START
|
|||
|
|
t_fit = df.loc[mask, "time_s"].values
|
|||
|
|
T_fit = df.loc[mask, "temp_obj IR [C]"].values
|
|||
|
|
|
|||
|
|
# Pesi espliciti: w=0 nella zona di transizione
|
|||
|
|
sigma = np.where(
|
|||
|
|
(t_fit >= W_ZERO_START) & (t_fit <= W_ZERO_END),
|
|||
|
|
1e10,
|
|||
|
|
1.0
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
# --- Modello doppio esponenziale ---
|
|||
|
|
def modello(t, A1, tau1, A2, tau2):
|
|||
|
|
return (T_INF
|
|||
|
|
+ A1 * np.exp(-(t - T1) / tau1)
|
|||
|
|
+ A2 * np.exp(-(t - T2) / tau2))
|
|||
|
|
|
|||
|
|
# Stime iniziali dai fit singoli
|
|||
|
|
p0 = [194.51, 13.17, 154.94, 17.12]
|
|||
|
|
bounds = ([0, 0.1, 0, 0.1], [np.inf, np.inf, np.inf, np.inf])
|
|||
|
|
|
|||
|
|
popt, pcov = curve_fit(
|
|||
|
|
modello, t_fit, T_fit,
|
|||
|
|
p0=p0,
|
|||
|
|
sigma=sigma,
|
|||
|
|
absolute_sigma=True,
|
|||
|
|
method="trf",
|
|||
|
|
bounds=bounds
|
|||
|
|
)
|
|||
|
|
A1_fit, tau1_fit, A2_fit, tau2_fit = popt
|
|||
|
|
perr = np.sqrt(np.diag(pcov))
|
|||
|
|
|
|||
|
|
# --- R² (solo punti con peso pieno) ---
|
|||
|
|
mask_w = (t_fit < W_ZERO_START) | (t_fit > W_ZERO_END)
|
|||
|
|
T_pred_w = modello(t_fit[mask_w], *popt)
|
|||
|
|
ss_res = np.sum((T_fit[mask_w] - T_pred_w) ** 2)
|
|||
|
|
ss_tot = np.sum((T_fit[mask_w] - T_fit[mask_w].mean()) ** 2)
|
|||
|
|
r2 = 1 - ss_res / ss_tot
|
|||
|
|
|
|||
|
|
print(f"A1 = {A1_fit:.4f} ± {perr[0]:.4f} °C")
|
|||
|
|
print(f"tau1 = {tau1_fit:.4f} ± {perr[1]:.4f} s")
|
|||
|
|
print(f"A2 = {A2_fit:.4f} ± {perr[2]:.4f} °C")
|
|||
|
|
print(f"tau2 = {tau2_fit:.4f} ± {perr[3]:.4f} s")
|
|||
|
|
print(f"R² = {r2:.6f} (punti con peso pieno)")
|
|||
|
|
|
|||
|
|
# --- Curve per il plot ---
|
|||
|
|
t_curve = np.linspace(T_START, df["time_s"].max(), 1000)
|
|||
|
|
T_tot = modello(t_curve, *popt)
|
|||
|
|
T_exp1 = T_INF + A1_fit * np.exp(-(t_curve - T1) / tau1_fit)
|
|||
|
|
T_exp2 = T_INF + A2_fit * np.exp(-(t_curve - T2) / tau2_fit)
|
|||
|
|
|
|||
|
|
# --- Plot ---
|
|||
|
|
fig, ax = plt.subplots(figsize=(12, 5))
|
|||
|
|
|
|||
|
|
ax.plot(t_fit, T_fit,
|
|||
|
|
color="steelblue", linewidth=0.8, label="Dati raw (temp_obj)")
|
|||
|
|
ax.axvspan(W_ZERO_START, W_ZERO_END,
|
|||
|
|
color="orange", alpha=0.25, label=f"Zona esclusa [{W_ZERO_START}–{W_ZERO_END} s]")
|
|||
|
|
ax.plot(t_curve, T_exp1,
|
|||
|
|
color="tomato", linewidth=1.2, linestyle=":",
|
|||
|
|
label=rf"$T_\infty + A_1 e^{{-(t-{T1})/\tau_1}}$ ($A_1$={A1_fit:.1f}, $\tau_1$={tau1_fit:.2f} s)")
|
|||
|
|
ax.plot(t_curve, T_exp2,
|
|||
|
|
color="seagreen", linewidth=1.2, linestyle=":",
|
|||
|
|
label=rf"$T_\infty + A_2 e^{{-(t-{T2})/\tau_2}}$ ($A_2$={A2_fit:.1f}, $\tau_2$={tau2_fit:.2f} s)")
|
|||
|
|
ax.plot(t_curve, T_tot,
|
|||
|
|
color="purple", linewidth=2, linestyle="--",
|
|||
|
|
label="Somma (fit totale)")
|
|||
|
|
|
|||
|
|
ax.set_xlabel("Tempo [s]")
|
|||
|
|
ax.set_ylabel("Temperatura [°C]")
|
|||
|
|
ax.set_title(f"Fit doppio esponenziale | R² = {r2:.4f}")
|
|||
|
|
ax.legend(fontsize=8)
|
|||
|
|
ax.grid(True, alpha=0.3)
|
|||
|
|
|
|||
|
|
plt.tight_layout()
|
|||
|
|
plt.savefig("fit_doppio_esponenziale.png", dpi=150, bbox_inches="tight")
|
|||
|
|
plt.show()
|