#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Parser inteligente de logs del WiT-CPV
(webpage_dynamictrack_autotune_log_psram.ino)
===============================================================================

Formato de lineas soportado:

    [AAAA-MM-DD HH:MM:SS] X:<x>,Y:<y>,C:<c>,V:<v>,P:<p>,S:<s>,er:<er>  -> muestra periodica (~2 s)
    [ts] FineTune: dC=.., dynE=.., e=.., near=0/1, Cmax=..             -> ciclo de fine-tune
    [ts] SunPos: Az=.., Alt=.., Az_corr=.., Alt_corr=.., posX=.., posY=.. -> prediccion astronomica
    [ts] Ejecutando comando: N                                         -> cambio de estado
    [ts] ErrMax=.. / MinErrMax=.. / CicleMax=.. / Az_deviation=..      -> parametros por web
    [ts] AT: rec small=.., abrupt=.., reset=.., peak_prox=..           -> auto-tune
    [ts] demas eventos (Timeinfo, Debug, Recibi, Moviendo, ...)

Calibracion (identica a metrica1.py), ADS122C04 de 23 bit, Vref = 2.048 V:

    G_POA [W/m2] = P_raw/8388607*2048/2      (piranometro, divisor /2)
    Voc   [V]    = V_raw/8388607*2.048*2     (divisor 10k/10k)
    Isc   [A]    = (C_raw/8388607*2048/20)/1e3/0.075   (AD8418 20x, shunt 75 mOhm)
    P_CPV [W]    = Voc * Isc * FF   con FF = 0.82

Geometria del soporte (trackSunByTime en el firmware):

    x = (90 - alt_corr)/80*53000 + 53000 + Az_deviation/90*y
    y = (az_corr - 180)/80*53000 + 53000
    => cada eje recorre 80 grados con 53000 micro-pasos (1 ustep = 80/53000 deg)

Salida: articulo/resultados_<tag>/{samples.csv, events.csv, metrics.json}
        articulo/figs/fig_{power,position,finetune,recovery}_<tag>.pdf

Uso:  py -3.12 parse_log.py 21092025.txt
"""

import json
import re
import sys
from pathlib import Path

import numpy as np
import pandas as pd

try:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.dates as mdates
    import matplotlib.pyplot as plt
    HAS_MPL = True
except Exception:
    HAS_MPL = False

# ----------------------------------------------------------------------------
# Calibracion (identica a metrica1.py)
# ----------------------------------------------------------------------------
ADC_MAX = 8388607.0        # 2^23 - 1
VREF_MV = 2048.0           # Vref interno del ADS122C04 (mV)
GAIN_SHUNT = 20.0          # ganancia del amplificador de instrumentacion
R_SHUNT = 0.075            # Ohm
FF = 0.82                  # fill factor caracterizado

# Geometria del soporte
STEPS_PER_80DEG = 53000.0  # micro-pasos por eje para 80 grados
DEG_PER_USTEP = 80.0 / STEPS_PER_80DEG

ROOT = Path(__file__).resolve().parent
ARTICULO = ROOT / "articulo"


def _trapz(y, x):
    """Integracion trapezoidal compatible con numpy >= 2 (np.trapezoid)."""
    if hasattr(np, "trapezoid"):
        return float(np.trapezoid(y, x=x))
    return float(np.trapz(y, x=x))


# ----------------------------------------------------------------------------
# Expresiones regulares de los eventos del log
# ----------------------------------------------------------------------------
RE_SAMPLE = re.compile(
    r"^X:(?P<x>-?\d+),Y:(?P<y>-?\d+),C:(?P<c>\d+),V:(?P<v>\d+),"
    r"P:(?P<p>\d+),S:(?P<s>\d+),er:(?P<er>\d+)\s*$"
)
RE_FINE = re.compile(
    r"FineTune: dC=(?P<dc>\d+), dynE=(?P<dynE>\d+), e=(?P<e>\d+), "
    r"near=(?P<near>\d), Cmax=(?P<cmax>\d+)"
)
RE_SUNPOS_NEW = re.compile(
    r"SunPos: Az=(?P<az>[\d.]+), Alt=(?P<alt>[\d.]+), "
    r"Az_corr=(?P<azc>[\d.]+), Alt_corr=(?P<altc>[\d.]+), "
    r"posX=(?P<px>\d+), posY=(?P<py>\d+)"
)
RE_SUNPOS_OLD = re.compile(
    r"SunPos: Az=(?P<az>[\d.]+), Zen=(?P<zen>[\d.]+) \(raw=[\d.]+\), "
    r"posX=(?P<px>\d+), posY=(?P<py>\d+)"
)
RE_CMD = re.compile(r"Ejecutando comando: (\d+)")
RE_PARAM = re.compile(r"^(ErrMax|MinErrMax|CicleMax|Az_deviation|TiltAngle)=(\S+)")
RE_AT_REC = re.compile(
    r"AT: rec small=(?P<small>\d+), abrupt=(?P<abrupt>\d+), "
    r"reset=(?P<reset>\d+), peak_prox=(?P<peak>[\d.]+)"
)
RE_AT_MEAN = re.compile(r"AT: meanI=([\d.]+)")


def parse_log(file_path):
    """Lee el log y devuelve (samples, events) como DataFrames."""
    samples = []
    events = []
    with open(file_path, "r", encoding="utf-8", errors="replace") as f:
        for line in f:
            line = line.rstrip("\n")
            # El firmware antiguo escribe cada linea con prefijo "LOG:" (antes
            # del timestamp). Retirarlo aqui (es un no-op en el formato nuevo)
            # permite que el mismo parser lea ambos formatos de log.
            if line.startswith("LOG:"):
                line = line[4:]
            m = re.match(r"^\[(\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2})\] (.*)$", line)
            if not m:
                continue
            ts, body = m.group(1), m.group(2)

            ms = RE_SAMPLE.match(body)
            if ms:
                samples.append(
                    {
                        "timestamp": ts,
                        "x": int(ms["x"]),
                        "y": int(ms["y"]),
                        "c_raw": int(ms["c"]),
                        "v_raw": int(ms["v"]),
                        "p_raw": int(ms["p"]),
                        "s_uptime": int(ms["s"]),
                        "er": int(ms["er"]),
                    }
                )
                continue

            mf = RE_FINE.search(body)
            if mf:
                events.append(
                    {"timestamp": ts, "kind": "fine",
                     "dC": int(mf["dc"]), "dynE": int(mf["dynE"]),
                     "e": int(mf["e"]), "near": int(mf["near"]),
                     "cmax": int(mf["cmax"])}
                )
                continue

            mp = RE_SUNPOS_NEW.search(body)
            if mp is None:
                mp = RE_SUNPOS_OLD.search(body)
            if mp:
                d = mp.groupdict()
                if "alt" in d and d["alt"] is not None:
                    alt, altc = float(d["alt"]), float(d["altc"])
                else:
                    zen = float(d["zen"])
                    alt = 90.0 - zen
                    altc = 90.0 - zen
                events.append(
                    {"timestamp": ts, "kind": "sunpos", "az": float(d["az"]),
                     "alt": alt, "azc": float(d.get("azc", d["az"])),
                     "altc": altc, "px": int(d["px"]), "py": int(d["py"])}
                )
                continue

            mc = RE_CMD.search(body)
            if mc:
                events.append({"timestamp": ts, "kind": "cmd", "cmd": int(mc.group(1))})
                continue

            mp = RE_PARAM.match(body)
            if mp:
                events.append({"timestamp": ts, "kind": "param",
                               "name": mp.group(1), "value": float(mp.group(2))})
                continue

            mat = RE_AT_REC.search(body)
            if mat:
                events.append({"timestamp": ts, "kind": "at_rec",
                                **{k: float(mat[k])
                                   for k in ("small", "abrupt", "reset", "peak")}})
                continue

            if body.startswith("AT:"):
                events.append({"timestamp": ts, "kind": "at_info", "text": body})
                continue

            # eventos menores (contexto)
            matched = False
            for prefix in ("Timeinfo:", "UT_hour:", "Debug:", "Recibi", "Moviendo"):
                if body.startswith(prefix):
                    events.append({"timestamp": ts, "kind": "info", "text": body})
                    matched = True
                    break
            if not matched:
                events.append({"timestamp": ts, "kind": "other", "text": body})

    df_s = pd.DataFrame(samples)
    if df_s.empty:
        raise ValueError("Sin muestras X,Y,C,V,P en el log.")
    df_s["timestamp"] = pd.to_datetime(df_s["timestamp"])
    df_e = pd.DataFrame(events)
    if not df_e.empty:
        df_e["timestamp"] = pd.to_datetime(df_e["timestamp"])
    return df_s, df_e



# ----------------------------------------------------------------------------
# Unidades fisicas + maquina de estados + referencia P_ref
# ----------------------------------------------------------------------------
CMD_NAMES = {2: "Track solar", 3: "Fine-tune", 4: "Homing",
              5: "Recentering", 6: "Set position", 10: "Driver status", 11: "Demo"}


def build_physical(df_s, df_e):
    """Añade unidades físicas, P_ref, estado/fase y posicion objetivo."""
    # 1) unidades fisicas (calibracion de metrica1.py)
    df_s["G_POA"] = df_s["p_raw"] / ADC_MAX * (VREF_MV / 2.0)          # W/m2
    df_s["Voc"] = df_s["v_raw"] / ADC_MAX * 2.048 * 2.0                # V
    df_s["Isc"] = (df_s["c_raw"] / ADC_MAX * (VREF_MV / GAIN_SHUNT)
                   / 1000.0) / R_SHUNT                                  # A
    df_s["P_CPV"] = df_s["Voc"] * df_s["Isc"] * FF                      # W

    # 2) maquina de estados a partir de los comandos recibidos
    status = 0
    phase = "Startup"
    phase_marks = []
    prev_ts = df_s["timestamp"].iloc[0]
    cmds = (df_e[df_e["kind"] == "cmd"] if not df_e.empty
            else pd.DataFrame(columns=["timestamp", "cmd"]))
    for ts_c, cmd in zip(cmds["timestamp"], cmds["cmd"]):
        phase_marks.append((prev_ts, ts_c, phase, status))
        if cmd in (2, 3):
            status = 3
            phase = "Solar tracking"
        elif cmd == 4:
            phase = "Homing"
        elif cmd == 5:
            phase = "Recentering"
        elif cmd == 11:
            phase = "Demo"
        else:
            phase = CMD_NAMES.get(cmd, f"Cmd {cmd}")
        prev_ts = ts_c
    phase_marks.append((prev_ts, df_s["timestamp"].iloc[-1], phase, status))

    # 3) asignar estado/fase a cada muestra
    ph = pd.DataFrame(phase_marks, columns=["t0", "t1", "phase", "status"])
    idx = ph["t0"].values.searchsorted(df_s["timestamp"].values, side="right") - 1
    idx = np.clip(idx, 0, len(ph) - 1)
    df_s["status"] = ph["status"].values[idx]
    df_s["phase"] = ph["phase"].values[idx]

    # 4) objetivo astronomico mas reciente (fwd-fill) + Ang_dev aplicado
    if not df_e.empty:
        sp = df_e[df_e["kind"] == "sunpos"].copy()
        if not sp.empty:
            tgt = pd.DataFrame({"timestamp": sp["timestamp"],
                                 "tgt_x": sp["px"], "tgt_y": sp["py"],
                                 "az": sp["az"], "alt": sp["alt"],
                                 "azc": sp["azc"], "altc": sp["altc"]})
            merged = pd.merge_asof(df_s[["timestamp"]], tgt,
                                    on="timestamp", direction="backward")
            for col in ("tgt_x", "tgt_y", "az", "alt", "azc", "altc"):
                df_s[col] = merged[col]
        azdev = 0.0
        par = df_e[df_e["kind"] == "param"]
        if not par.empty:
            azp = par[par["name"] == "Az_deviation"]
            if not azp.empty:
                azser = pd.Series(azp["value"].values, index=azp["timestamp"])
                ff = azser.reindex(df_s["timestamp"], method="ffill").dropna()
                azdev = float(ff.iloc[0]) if len(ff) else 0.0
        # angulos reconstruidos invirtiendo el mapeo del firmware
        x_base = df_s["x"] - azdev / 90.0 * df_s["y"]
        df_s["alt_est"] = 90.0 - (x_base - STEPS_PER_80DEG) / STEPS_PER_80DEG * 80.0
        df_s["az_est"] = 180.0 + (df_s["y"] - STEPS_PER_80DEG) / STEPS_PER_80DEG * 80.0
        df_s["azdev_used"] = azdev
    else:
        for col in ("tgt_x", "tgt_y", "az", "alt", "azc", "altc", "azdev_used"):
            df_s[col] = np.nan
        df_s["azdev_used"] = 0.0

    # 5) ciclos de fine-tune alineados por timestamp (merge asof)
    if not df_e.empty:
        ft = df_e[df_e["kind"] == "fine"].copy()
        if not ft.empty:
            df_s = pd.merge_asof(df_s, ft[["timestamp", "e", "dC", "dynE",
                                            "near", "cmax"]].sort_values("timestamp"),
                                  on="timestamp", direction="backward")
    else:
        for col in ("e", "dC", "dynE", "near", "cmax"):
            df_s[col] = np.nan

    # 6) referencia optima: P_ref(t) = R_opt,max * G(t) * Voc(t) * FF
    #    R_opt,max se estima bajo cielo directo (G >= 600 W/m2) para que los
    #    transitorios nubosos (respuesta angular del piranometro vs optica
    #    concentradora) no inflen la responsividad de referencia.
    lit = df_s[df_s["G_POA"] > 600.0]
    R_opt = float((lit["Isc"] / lit["G_POA"]).max()) if len(lit) else np.nan
    df_s["R_opt"] = R_opt
    df_s["Isc_ref"] = R_opt * df_s["G_POA"]
    df_s["P_ref"] = R_opt * df_s["G_POA"] * df_s["Voc"] * FF

    # 7) actividad mecanica e intervalos temporales
    df_s["dx"] = df_s["x"].diff().abs()
    df_s["dy"] = df_s["y"].diff().abs()
    df_s["dt_s"] = df_s["timestamp"].diff().dt.total_seconds().fillna(0.0)

    return df_s, ph



# ----------------------------------------------------------------------------
# Metricas (seccion 4 del articulo)
# ----------------------------------------------------------------------------
def window_metrics(df, label):
    """Metricas de una ventana del log (df ya enriquecido con P_ref)."""
    w = df[df["G_POA"] > 50.0] if len(df) else df
    if w.empty:
        return {"label": label, "n_samples": 0}
    # segundos desde el inicio de la ventana (independiente de la resolucion
    # de datetime de pandas 3.x)
    ts = w["timestamp"].values
    x = (ts - ts[0]) / np.timedelta64(1, "s")
    E = _trapz(w["P_CPV"].values, x) / 3600.0            # Wh
    Pref = _trapz(w["P_ref"].values, x)                  # W.s
    H = _trapz(w["G_POA"].values, x) / 3600.0 / 1000.0   # kWh/m2
    eta = 100.0 * E * 3600.0 / Pref if Pref > 0 else np.nan
    S_el, S_az = int(w["dx"].sum()), int(w["dy"].sum())
    return {
        "label": label,
        "n_samples": int(len(w)),
        "t_inicio": str(w["timestamp"].iloc[0]),
        "t_fin": str(w["timestamp"].iloc[-1]),
        "E_Wh": round(E, 3),
        "H_kWh_m2": round(H, 4),
        "E_norm_Wh_per_kWhm2": round(E / H, 3) if H > 0 else None,
        "eta_tr_pct": round(eta, 2) if eta == eta else None,
        "MAE_I_mA": round(1000.0 * np.mean(np.abs(w["Isc_ref"] - w["Isc"])), 2),
        "MAE_P_mW": round(1000.0 * np.mean(np.abs(w["P_ref"] - w["P_CPV"])), 2),
        "P_CPV_mean_mW": round(1000.0 * np.mean(w["P_CPV"]), 1),
        "P_ref_mean_mW": round(1000.0 * np.mean(w["P_ref"]), 1),
        "S_el_usteps": S_el,
        "S_az_usteps": S_az,
        "S_tot_usteps": S_el + S_az,
        "S_tot_deg": round((S_el + S_az) * DEG_PER_USTEP, 1),
    }


def recovery_t95(df, df_e, horizon_min=20):
    """t95 tras cada comando de re-apuntado/re-centrado.

    Se define la recuperacion como el primer instante, en una ventana de
    `horizon_min` minutos despues del trigger, en que P_CPV alcanza de nuevo
    el 95% de P_ref. Si no se alcanza dentro de la ventana se reporta None
    (y, como dato informativo, el instante en que si se alcanza).
    """
    out = []
    cmds = df_e[df_e["kind"] == "cmd"] if not df_e.empty else pd.DataFrame()
    if cmds.empty:
        return pd.DataFrame(out)
    for _, row in cmds.iterrows():
        n, t0 = int(row["cmd"]), row["timestamp"]
        if n not in (2, 5):
            continue
        w = df[(df["timestamp"] >= t0) &
               (df["timestamp"] <= t0 + pd.Timedelta(minutes=horizon_min))].reset_index(drop=True)
        if w.empty:
            continue
        thr = 0.95 * w["P_ref"].values
        P = w["P_CPV"].values
        below = np.where(P < thr)[0]
        if below.size == 0:
            out.append({"trigger": CMD_NAMES.get(n, n), "t_trigger": str(t0),
                        "t95_s": 0.0, "dropped": False, "P_min_mW": None})
            continue
        i_drop = int(below[0])
        P_min = float(P[i_drop:].min())
        t95 = None
        above = np.where(P[i_drop:] >= thr)[0]
        if above.size:
            t95 = (w["timestamp"].iloc[i_drop + int(above[0])] - t0).total_seconds()
        # cruce posterior (fuera de la ventana), solo informativo
        t_late = None
        rest = df[df["timestamp"] > t0 + pd.Timedelta(minutes=horizon_min)]
        if t95 is None and not rest.empty:
            above_late = np.where(rest["P_CPV"].values >= 0.95 * rest["P_ref"].values)[0]
            if above_late.size:
                t_late = (rest["timestamp"].iloc[int(above_late[0])] - t0).total_seconds()
        out.append({"trigger": CMD_NAMES.get(n, n), "t_trigger": str(t0),
                    "t95_s": round(t95, 1) if t95 is not None else None,
                    "t95_late_s": round(t_late, 1) if t_late is not None else None,
                    "dropped": True, "P_min_mW": round(P_min * 1000, 1)})
    return pd.DataFrame(out)


def cloud_events(df):
    """Transitorios nubosos: G < 0.85 * mediana movil 10 min (ventanas >= 120 s)."""
    if df.empty:
        return []
    d = df.set_index("timestamp").sort_index()
    g = d["G_POA"]
    if g.max() < 300:
        return []
    med = g.rolling("10min", center=True, min_periods=10).median()
    mask = (g < 0.85 * med) & (g > 30)
    events = []
    t = d.index
    blocks, start = [], None
    for i in range(len(mask)):
        if mask.iloc[i] and start is None:
            start = i
        elif not mask.iloc[i] and start is not None:
            blocks.append((start, i - 1))
            start = None
    if start is not None:
        blocks.append((start, len(mask) - 1))
    for i0, i1 in blocks:
        dur = (t[i1] - t[i0]).total_seconds()
        if dur < 120:
            continue
        g_seg = g.iloc[i0:i1 + 1]
        base = float(g.iloc[max(0, i0 - 5):i0].mean())
        g_min = float(g_seg.min())
        rec_idx = np.where(g.iloc[i1:].values >= 0.95 * base)[0]
        t_rec = (t[i1 + int(rec_idx[0])] - t[i0]).total_seconds() if rec_idx.size else None
        events.append({
            "t_start": str(t[i0]), "dur_s": int(dur),
            "G_base_Wm2": round(base, 0), "G_min_Wm2": round(g_min, 0),
            "depth_pct": round(100.0 * (1 - g_min / base), 1) if base > 0 else None,
            "t95_G_s": round(t_rec, 1) if t_rec is not None else None,
        })
    return events


def position_errors(df, df_e):
    """Error de posicion vs objetivo astronomico tras cada reticle.

    'inicial': primera muestra despues de la orden (fin del recorrido).
    'convergido': mediana del offset durante los ultimos 5 min de los
    primeros 10 min tras el reticle (periodo de convergencia del fine-tune;
    despues, el objetivo geometrico queda obsoleto porque el sol sigue
    moviendose y el controlador persigue el pico electrico, no el target).
    """
    out = []
    if df_e.empty:
        return pd.DataFrame(out)
    sp = df_e[df_e["kind"] == "sunpos"]
    for i, (_, row) in enumerate(sp.iterrows()):
        t0, tx, ty = row["timestamp"], row["px"], row["py"]
        seg = df[(df["timestamp"] >= t0) & (df["status"] == 3)]
        if seg.empty:
            continue
        first = seg.iloc[0]
        conv = seg[(seg["timestamp"] >= t0 + pd.Timedelta(minutes=5)) &
                   (seg["timestamp"] <= t0 + pd.Timedelta(minutes=10))]
        if conv.empty:
            conv = seg.tail(50)
        ex0, ey0 = abs(first["x"] - tx) * DEG_PER_USTEP, abs(first["y"] - ty) * DEG_PER_USTEP
        exs = float(np.median([abs(r["x"] - tx) * DEG_PER_USTEP for _, r in conv.iterrows()]))
        eys = float(np.median([abs(r["y"] - ty) * DEG_PER_USTEP for _, r in conv.iterrows()]))
        out.append({
            "t_reticle": str(t0), "az_corr": row["azc"], "alt_corr": row["altc"],
            "errX0_deg": round(ex0, 2), "errY0_deg": round(ey0, 2),
            "errX0_usteps": int(abs(first["x"] - tx)), "errY0_usteps": int(abs(first["y"] - ty)),
            "errX_ss_deg": round(exs, 3), "errY_ss_deg": round(eys, 3),
        })
    return pd.DataFrame(out)



# ----------------------------------------------------------------------------
# Figuras (PDF vectoriales para el articulo)
# ----------------------------------------------------------------------------
CMD_COLORS = {2: "tab:green", 3: "tab:purple", 4: "tab:blue", 5: "tab:brown"}


def make_figures(df, df_e, outdir, tag):
    if not HAS_MPL:
        print("[aviso] matplotlib no disponible; se omiten figuras.")
        return
    figs_dir = ARTICULO / "figs"
    figs_dir.mkdir(parents=True, exist_ok=True)
    t = df["timestamp"]
    lit_mask = (df["G_POA"] > 50.0).values
    cmds = df_e[df_e["kind"] == "cmd"] if not df_e.empty else pd.DataFrame()

    # --- Fig 1: potencia generada vs referencia + irradiancia
    fig, ax1 = plt.subplots(figsize=(7.2, 3.4), dpi=150)
    ax1.plot(t, df["P_CPV"] * 1000, color="tab:red", lw=1.0, label="$P_{CPV}$")
    ax1.plot(t, df["P_ref"] * 1000, color="k", ls="--", lw=1.0, label="$P_{ref}$")
    ax1.fill_between(t, df["P_CPV"] * 1000, df["P_ref"] * 1000,
                     where=lit_mask, color="tab:red", alpha=0.12,
                     label="Loss vs. $P_{ref}$")
    ax1.set_ylabel("Power (mW)", color="tab:red")
    ax1.tick_params(axis="y", labelcolor="tab:red")
    ax2 = ax1.twinx()
    ax2.plot(t, df["G_POA"], color="tab:orange", lw=0.8, alpha=0.85, label="$G_{POA}$")
    ax2.set_ylabel("Irradiance (W/m$^2$)", color="tab:orange")
    ax2.tick_params(axis="y", labelcolor="tab:orange")
    labeled = False
    for _, r in cmds.iterrows():
        if int(r["cmd"]) != 2:
            continue  # solo se anota "Track solar": los demas comandos
        # estan todos dentro de ~75 s entre si y sus etiquetas se solaparian
        ax1.axvline(r["timestamp"], color=CMD_COLORS.get(r["cmd"], "gray"),
                    ls=":", lw=1.0, alpha=0.8)
        if not labeled:  # etiqueta unica: ambos Track-solar estan a < 2 min
            ax1.text(r["timestamp"], ax1.get_ylim()[1] * 0.95,
                     f"cmd {int(r['cmd'])}\n({CMD_NAMES.get(int(r['cmd']), '')})",
                     fontsize=6.5, color=CMD_COLORS.get(r["cmd"], "gray"))
            labeled = True
    h1, l1 = ax1.get_legend_handles_labels()
    h2, l2 = ax2.get_legend_handles_labels()
    ax1.legend(h1 + h2, l1 + l2, loc="upper right", fontsize=7, framealpha=0.6)
    ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    ax1.set_xlabel("Local time (21 Sep 2025, UTC-6)")
    ax1.set_title(f"WiT-CPV field test - session {tag}", fontsize=9)
    fig.tight_layout()
    fig.savefig(figs_dir / f"fig_power_{tag}.pdf")
    plt.close(fig)

    # --- Fig 2: posicion real vs objetivo astronomico
    fig, axes = plt.subplots(2, 1, figsize=(7.2, 4.2), dpi=150, sharex=True)
    for ax, col, name, color in ((axes[0], "x", "$X$ (zenithal, $\\mu$steps)", "tab:blue"),
                                  (axes[1], "y", "$Y$ (azimuthal, $\\mu$steps)", "tab:olive")):
        ax.plot(t, df[col], color=color, lw=0.9, label="Actual position")
        tc = df["tgt_x" if col == "x" else "tgt_y"]
        ax.step(t, tc, where="post", color="k", ls="--", lw=1.0,
                label="Predicted sun position")
        ax.set_ylabel(name)
        if not cmds.empty:
            for _, r in cmds.iterrows():
                ax.axvline(r["timestamp"], color=CMD_COLORS.get(r["cmd"], "gray"),
                           ls=":", lw=0.8, alpha=0.7)
        ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
        ax.tick_params(labelsize=7)
    axes[1].set_xlabel("Local time (21 Sep 2025, UTC-6)")
    fig.tight_layout()
    fig.savefig(figs_dir / f"fig_position_{tag}.pdf")
    plt.close(fig)



    # --- Fig 3: fine-tune adaptativo (e_dyn, |dC|)
    fig, ax1 = plt.subplots(figsize=(7.2, 3.2), dpi=150)
    ax1.plot(t, df["er"], color="tab:purple", lw=0.7, alpha=0.9,
              label="Dynamic threshold $e_{dyn}$")
    ax2 = ax1.twinx()
    ax2.semilogy(t, df["dC"].clip(lower=1), color="tab:red", lw=0.5, alpha=0.7,
                 label="$|\\Delta I_{sc}|$")
    ax2.set_ylabel("$|\\Delta I_{sc}|$ (ADC counts, log)", color="tab:red")
    ax2.tick_params(axis="y", labelcolor="tab:red")
    ax1.set_ylabel("$e_{dyn}$ ($\\mu$steps)")
    at = df_e[df_e["kind"] == "at_rec"] if not df_e.empty else pd.DataFrame()
    for _, r in at.iterrows():
        ax1.axvline(r["timestamp"], color="k", ls=":", lw=1.0)
        ax1.text(r["timestamp"], ax1.get_ylim()[1] * 0.9, "auto-tune\nrecommendation",
                 fontsize=6.5, ha="left")
    h1, l1 = ax1.get_legend_handles_labels()
    h2, l2 = ax2.get_legend_handles_labels()
    ax1.legend(h1 + h2, l1 + l2, loc="upper right", fontsize=7)
    ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    ax1.set_xlabel("Local time (21 Sep 2025, UTC-6)")
    fig.tight_layout()
    fig.savefig(figs_dir / f"fig_finetune_{tag}.pdf")
    plt.close(fig)

    # --- Fig 4: recuperacion t95 alrededor de los triggers
    rec = [r for _, r in cmds.iterrows() if int(r["cmd"]) in (2, 5)]
    if rec:
        fig, axes = plt.subplots(1, min(2, len(rec)), figsize=(7.2, 2.9), dpi=150)
        axes = np.atleast_1d(axes)
        for ax, r in zip(axes, rec[:2]):
            t0 = r["timestamp"]
            w = df[(df["timestamp"] >= t0 - pd.Timedelta(minutes=1)) &
                   (df["timestamp"] <= t0 + pd.Timedelta(minutes=20))]
            ax.plot(w["timestamp"], w["P_CPV"] * 1000, color="tab:red", lw=1.0)
            ax.plot(w["timestamp"], w["P_ref"] * 1000, color="k", ls="--", lw=1.0)
            ax.plot(w["timestamp"], 0.95 * w["P_ref"] * 1000, color="tab:green",
                    ls=":", lw=1.0, label="0.95 $P_{ref}$")
            ax.axvline(t0, color="tab:blue", ls=":", lw=1.0)
            thr = 0.95 * w["P_ref"].values
            P = w["P_CPV"].values
            below = np.where(P < thr)[0]
            if below.size:
                i = int(below[0])
                above = np.where(P[i:] >= thr)[0]
                if above.size:
                    t95s = (w["timestamp"].iloc[i + int(above[0])] -
                            w["timestamp"].iloc[i]).total_seconds()
                    ax.annotate(f"$t_{{95}}$ = {t95s:.0f} s",
                                xy=(w["timestamp"].iloc[i + int(above[0])],
                                    P[i + int(above[0])] * 1000),
                                fontsize=8, color="tab:green")
            ax.set_title(f"cmd {int(r['cmd'])} ({CMD_NAMES.get(int(r['cmd']), '')}) "
                         f"@ {str(r['timestamp'])[11:16]}", fontsize=8)
            ax.set_ylabel("Power (mW)")
            ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))
            ax.tick_params(labelsize=7)
            ax.legend(fontsize=7)
        fig.tight_layout()
        fig.savefig(figs_dir / f"fig_recovery_{tag}.pdf")
        plt.close(fig)
    print(f"[figs] guardadas en {figs_dir}")



# ----------------------------------------------------------------------------
# Main
# ----------------------------------------------------------------------------
def main():
    args = sys.argv[1:]
    log_file = args[0] if args else "21092025.txt"
    tag = Path(log_file).stem
    path = ROOT / log_file
    outdir = ARTICULO / f"resultados_{tag}"
    outdir.mkdir(parents=True, exist_ok=True)

    print(f"[1/5] Parseando {path} ...")
    df_s, df_e = parse_log(path)
    print(f"      muestras={len(df_s)}  eventos={len(df_e)}  "
          f"({df_s['timestamp'].iloc[0]} -> {df_s['timestamp'].iloc[-1]})")

    print("[2/5] Unidades fisicas + estados + P_ref ...")
    df, ph = build_physical(df_s, df_e)
    lit = df[df["G_POA"] > 200.0]
    R_opt = df["R_opt"].iloc[0]
    print(f"      R_opt,max = {1000 * R_opt:.3f} mA/(W/m2)   "
          f"(ventana iluminada: {lit['timestamp'].iloc[0]} -> {lit['timestamp'].iloc[-1]})")

    print("[3/5] Metricas ...")
    phases = ph[ph["status"] == 3]
    metrics = {"tag": tag,
               "log": path.name,
               "session_start": str(df["timestamp"].iloc[0]),
               "session_end": str(df["timestamp"].iloc[-1]),
               "R_opt_max_mA_per_Wm2": round(1000 * R_opt, 3),
               "FF_used": FF}
    metrics["full_session"] = window_metrics(df, "Full session (G > 50 W/m2)")
    metrics["daylight"] = window_metrics(lit, "Daylight (G > 200 W/m2)")
    metrics["tracking_active"] = window_metrics(df[df["status"] == 3],
                                                "Tracking active (status=3)")
    clear = df[df["G_POA"] > 600.0]
    if len(clear):
        m_clear = window_metrics(clear, "Clear-sky (G > 600 W/m2)")
        m_clear["frac_P_above_095Pref"] = round(
            float((clear["P_CPV"] >= 0.95 * clear["P_ref"]).mean()), 3)
        m_clear["ratio_P_Pref_median"] = round(
            float((clear["P_CPV"] / clear["P_ref"]).median()), 3)
        metrics["clear_sky"] = m_clear
    else:
        metrics["clear_sky"] = None
    metrics["phases"] = [window_metrics(df[df["phase"] == p], f"Phase: {p}")
                         for p in phases["phase"].drop_duplicates()]

    rec = recovery_t95(df, df_e)
    metrics["t95_recovery"] = rec.to_dict("records") if not rec.empty else []
    metrics["cloud_events"] = cloud_events(df)
    pe = position_errors(df, df_e)
    metrics["position_errors"] = pe.to_dict("records") if not pe.empty else []
    er_seq = df["er"].dropna().astype(int)
    changes = er_seq[er_seq != er_seq.shift(1)]
    changes = changes.iloc[1:] if changes.index[0] == df.index[0] else changes
    metrics["er_progression"] = [
        {"t": str(ts), "e": int(ev)}
        for ts, ev in zip(df["timestamp"].values[changes.index], changes.values)]
    if not df_e.empty:
        at = df_e[df_e["kind"] == "at_rec"]
        cols = [c for c in ("timestamp", "small", "abrupt", "reset", "peak")
                if c in at.columns]
        metrics["at_recommendations"] = at[cols].to_dict("records")
    else:
        metrics["at_recommendations"] = []

    print("[4/5] Guardando salida ...")
    keep = ["timestamp", "x", "y", "c_raw", "v_raw", "p_raw", "s_uptime", "er",
            "G_POA", "Voc", "Isc", "P_CPV", "Isc_ref", "P_ref", "status", "phase",
            "az_est", "alt_est", "tgt_x", "tgt_y", "az", "alt", "azc", "altc",
            "azdev_used", "e", "dC", "dynE", "near", "cmax"]
    df[keep].to_csv(outdir / "samples.csv", index=False)
    df_e.to_csv(outdir / "events.csv", index=False, na_rep="")
    with open(outdir / "metrics.json", "w", encoding="utf-8") as f:
        json.dump(metrics, f, indent=2, ensure_ascii=False, default=str)

    print("[5/5] Figuras ...")
    make_figures(df, df_e, outdir, tag)

    # -------- resumen en consola --------
    print("\n" + "=" * 68)
    print(f"   RESUMEN DE METRICAS  (log {path.name})")
    print("=" * 68)
    for key in ("full_session", "daylight", "tracking_active", "clear_sky"):
        m = metrics[key]
        if m is None:
            continue
        extra = (f" | %P>=0.95Pref = {m['frac_P_above_095Pref']:.0%} | "
                 f"med(P/Pref) = {m['ratio_P_Pref_median']}"
                 if "frac_P_above_095Pref" in m else "")
        print(f"\n--- {m['label']}  ({m['t_inicio'][11:19]} -> {m['t_fin'][11:19]}, "
              f"n={m['n_samples']})")
        print(f"    E = {m['E_Wh']:.2f} Wh | H = {m['H_kWh_m2']:.4f} kWh/m2 | "
              f"E/H = {m['E_norm_Wh_per_kWhm2']} Wh/(kWh/m2)")
        print(f"    eta_tr = {m['eta_tr_pct']} % | MAE_I = {m['MAE_I_mA']} mA | "
              f"MAE_P = {m['MAE_P_mW']} mW{extra}")
        print(f"    S_el = {m['S_el_usteps']:,} | S_az = {m['S_az_usteps']:,} | "
              f"S_tot = {m['S_tot_usteps']:,} usteps ({m['S_tot_deg']} deg)")
    for m in metrics["phases"]:
        print(f"\n    [{m['label']}] {m['t_inicio'][11:16]} -> {m['t_fin'][11:16]} | "
              f"eta={m['eta_tr_pct']}% E={m['E_Wh']}Wh H={m['H_kWh_m2']}kWh/m2 "
              f"S={m['S_tot_usteps']:,}")
    if metrics["t95_recovery"]:
        print("\n--- t95 (recuperacion a 95% de P_ref tras trigger, horizonte 20 min)")
        for r in metrics["t95_recovery"]:
            extra = (f" (P_min = {r['P_min_mW']} mW)"
                     if r.get("P_min_mW") is not None else "")
            late = (f" [primer cruce posterior: {r['t95_late_s']} s despues]"
                    if r.get("t95_late_s") else "")
            print(f"    {r['trigger']:>12} @ {r['t_trigger'][11:19]} -> t95 = "
                  f"{r['t95_s']} s{extra}{late}")
    if metrics["cloud_events"]:
        print(f"\n--- {len(metrics['cloud_events'])} transitorios nubosos detectados")
        for c in metrics["cloud_events"][:8]:
            print(f"    {c['t_start'][11:19]}  G: {c['G_base_Wm2']:.0f} -> "
                  f"{c['G_min_Wm2']:.0f} W/m2 ({c['depth_pct']}%)  "
                  f"dur={c['dur_s']}s  t95_G={c['t95_G_s']}s")
    if metrics["position_errors"]:
        print("\n--- Error de posicion vs objetivo astronomico")
        for p in metrics["position_errors"]:
            print(f"    {p['t_reticle'][11:19]}  inicial: {p['errX0_deg']}/"
                  f"{p['errY0_deg']} deg ({p['errX0_usteps']}/{p['errY0_usteps']} "
                  f"usteps) | regimen: {p['errX_ss_deg']}/{p['errY_ss_deg']} deg")
    if metrics.get("er_progression"):
        print("\n--- Progresion del umbral adaptativo e_dyn (cambios)")
        print("    " + " -> ".join(f"{r['t'][11:19]}:{r['e']}"
                                    for r in metrics["er_progression"]))
    if metrics.get("at_recommendations"):
        for a in metrics["at_recommendations"]:
            print(f"\n--- Auto-tune @ {str(a['timestamp'])[11:19]}: "
                  f"small={a['small']:.0f}, abrupt={a['abrupt']:.0f}, "
                  f"reset={a['reset']:.0f}, peak_prox={a['peak']}")
    print(f"\n[ok] salida en {outdir}\n")


if __name__ == "__main__":
    main()