import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from pathlib import Path

# --- Dati ---
ROOT = Path(__file__).parent.parent
df = pd.read_csv(ROOT / "data.csv")
df["time_s"] = df["time since start [ms]"] / 1000.0

T_INF    = 22.99   # temperatura ambiente [°C]
T0       = 115.0   # inizio intervallo [s]
T_END    = 133.7   # fine intervallo [s] — ultimo campione nel dataset
W0_START = 115.9   # inizio zona esclusa [s]  (None = nessuna esclusione)
W0_END   = 117.2   # fine zona esclusa [s]

mask  = (df["time_s"] >= T0) & (df["time_s"] <= T_END)
t_fit = df.loc[mask, "time_s"].values
T_fit = df.loc[mask, "temp_obj IR [C]"].values

sigma = np.where(
    (t_fit >= W0_START) & (t_fit <= W0_END),
    1e10,
    1.0
)

# --- Modello: esponenziale di Newton ---
def modello(t, A, tau):
    return T_INF + A * np.exp(-(t - T0) / tau)

A0   = T_fit[0] - T_INF
tau0 = 5.0

popt, pcov = curve_fit(
    modello, t_fit, T_fit,
    p0=[A0, tau0],
    sigma=sigma, absolute_sigma=True,
    method="trf",
    bounds=([-np.inf, 0.01], [np.inf, np.inf])
)
A_fit, tau_fit = popt
perr = np.sqrt(np.diag(pcov))

# --- R² (solo punti con peso pieno) ---
mask_clean = (t_fit < W0_START) | (t_fit > W0_END)
T_pred = modello(t_fit[mask_clean], *popt)
ss_res = np.sum((T_fit[mask_clean] - T_pred) ** 2)
ss_tot = np.sum((T_fit[mask_clean] - T_fit[mask_clean].mean()) ** 2)
r2 = 1 - ss_res / ss_tot

print(f"A   = {A_fit:.4f} +/- {perr[0]:.4f} °C")
print(f"tau = {tau_fit:.4f} +/- {perr[1]:.4f} s")
print(f"R²  = {r2:.6f}")
print(f"\nEquazione: T(t) = {T_INF} + {A_fit:.4f} * exp(-(t - {T0}) / {tau_fit:.4f})")

# --- Curva continua per il plot ---
t_curve = np.linspace(T0, T_END, 500)
T_curve = modello(t_curve, *popt)

# --- Plot ---
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(t_fit, T_fit,
           color="steelblue", s=20, zorder=3, label="Dati raw (temp_obj)")
ax.axvspan(W0_START, W0_END, color="orange", alpha=0.25,
           label=f"Zona esclusa [{W0_START}–{W0_END} s]")
ax.plot(t_curve, T_curve,
        color="tomato", linewidth=2,
        label=f"Fit: {T_INF} + {A_fit:.2f}·exp(-(t-{T0})/{tau_fit:.3f})")
ax.set_xlabel("Tempo [s]")
ax.set_ylabel("Temperatura [°C]")
ax.set_title(f"Fit esponenziale  [{T0}–{T_END} s]  |  R² = {r2:.4f}")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(Path(__file__).parent / "induttore3.png", dpi=150, bbox_inches="tight")
plt.show()

# --- Report ---
report = Path(__file__).parent / "report3.md"
with open(report, "w", encoding="utf-8") as f:
    f.write("# Report Induttore 3\n\n")
    f.write("## Intervallo di campionamento\n\n")
    f.write(f"| Parametro | Valore |\n|---|---|\n")
    f.write(f"| Inizio | {T0} s |\n")
    f.write(f"| Fine | {T_END} s |\n")
    f.write(f"| Numero di punti | {len(t_fit)} |\n")
    f.write(f"| Zona esclusa | {W0_START}–{W0_END} s |\n\n")
    f.write("## Equazione trovata\n\n")
    f.write(f"$$T(t) = {T_INF} + {A_fit:.4f} \\cdot e^{{-(t - {T0}) / {tau_fit:.4f}}}$$\n\n")
    f.write("| Parametro | Valore | Incertezza (1σ) |\n|---|---|---|\n")
    f.write(f"| A | {A_fit:.4f} °C | +/- {perr[0]:.4f} °C |\n")
    f.write(f"| τ | {tau_fit:.4f} s | +/- {perr[1]:.4f} s |\n\n")
    f.write("## Qualità del fit\n\n")
    f.write(f"| R² | {r2:.6f} |\n|---|---|\n")
print("report3.md generato.")