"""Reproducible read-and-report analysis of the 4 November 2026 field log."""

from __future__ import annotations

import hashlib
import json
import re
import sys
from datetime import datetime
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd


BASE = Path(__file__).resolve().parent
SOURCE = BASE / "logs_nuevos" / "log_2026-11-04.txt"
SUNCALC_SOURCE = BASE / "logs_nuevos" / "suncalc_log_2026-11-04.pdf"
OUT = BASE / "figures" / "2026-11-04"
RESULT_TEX = BASE / "articulo" / "results_2026_11_04.tex"
METRICS_TEX = BASE / "articulo" / "metrics_2026_11_04.tex"

ADC_MAX = 8_388_607
VREF_V = 2.048
CURRENT_GAIN = 20.0
SHUNT_OHM = 0.130
PYRANOMETER_MV_PER_W_M2 = 0.2
APERTURE_M2 = 0.002497
FILL_FACTOR = 0.82
MAX_GAP_S = 15.0
MIN_IRRADIANCE_W_M2 = 50.0
SAMPLE_RE = re.compile(
    r"LOG:\[(?P<timestamp>[^\]]+)\] X:(?P<x>-?\d+),Y:(?P<y>-?\d+),"
    r"C:(?P<c>-?\d+),V:(?P<v>-?\d+),P:(?P<p>-?\d+),"
    r"S:(?P<uptime>\d+),er:(?P<error>\d+)"
)
TEL_RE = re.compile(
    r"LOG:\[(?P<timestamp>[^\]]+)\] TEL: t=(?P<uptime>\d+)s, "
    r"steps=(?P<steps>\d+)\(\+?(?P<delta>\d+)\)"
)
EVENT_RE = re.compile(r"LOG:\[(?P<timestamp>[^\]]+)\] (?P<message>.*)")

# Exact 15-minute values transcribed from the supplied SunCalc PDF.
SUNCALC_ROWS = [
    ("07:00",106.87,.82),("07:15",108.31,3.91),("07:30",109.82,7.12),("07:45",111.41,10.34),
    ("08:00",113.08,13.54),("08:15",114.86,16.70),("08:30",116.76,19.83),("08:45",118.80,22.90),
    ("09:00",120.98,25.91),("09:15",123.34,28.85),("09:30",125.90,31.71),("09:45",128.67,34.47),
    ("10:00",131.69,37.13),("10:15",134.98,39.66),("10:30",138.57,42.04),("10:45",142.49,44.24),
    ("11:00",146.75,46.26),("11:15",151.37,48.04),("11:30",156.34,49.57),("11:45",161.64,50.82),
    ("12:00",167.22,51.75),("12:15",173.01,52.35),("12:30",178.93,52.59),("12:45",184.86,52.47),
    ("13:00",190.71,51.99),("13:15",196.37,51.18),("13:30",201.77,50.04),("13:45",206.86,48.60),
    ("14:00",211.60,46.90),("14:15",215.98,44.96),("14:30",220.01,42.81),("14:45",223.71,40.49),
    ("15:00",227.10,38.01),("15:15",230.21,35.40),("15:30",233.06,32.67),("15:45",235.69,29.84),
    ("16:00",238.11,26.93),("16:15",240.35,23.94),("16:30",242.43,20.89),("16:45",244.37,17.78),
    ("17:00",246.18,14.63),("17:15",247.89,11.44),("17:30",249.50,8.23),("17:45",251.03,5.01),
    ("18:00",252.49,1.85),("18:15",253.88,-1.05),
]


def signed24(value: int) -> int:
    value &= 0xFFFFFF
    return value - 0x1000000 if value & 0x800000 else value


def parse_log(path: Path = SOURCE):
    samples, telemetry, events = [], [], []
    raw = path.read_bytes()
    for line_number, line in enumerate(raw.decode("utf-8-sig").splitlines(), 1):
        match = SAMPLE_RE.search(line)
        if match:
            d = match.groupdict()
            samples.append({
                "timestamp": pd.Timestamp(d["timestamp"]),
                "x_steps": int(d["x"]), "y_steps": int(d["y"]),
                "c_raw": int(d["c"]), "v_raw": int(d["v"]), "p_raw": int(d["p"]),
                "uptime_s": int(d["uptime"]), "error": int(d["error"]),
                "source_line": line_number,
            })
            continue
        match = TEL_RE.search(line)
        if match:
            d = match.groupdict()
            telemetry.append({"timestamp": pd.Timestamp(d["timestamp"]),
                              "uptime_s": int(d["uptime"]), "steps": int(d["steps"]),
                              "delta": int(d["delta"]), "source_line": line_number})
            continue
        match = EVENT_RE.search(line)
        if match:
            events.append({"timestamp": pd.Timestamp(match["timestamp"]),
                           "message": match["message"], "source_line": line_number})
    frame = pd.DataFrame(samples)
    if frame.empty or not frame.timestamp.is_monotonic_increasing:
        raise ValueError("No electrical samples or non-monotonic source timestamps")
    frame.attrs["firmware_versions"] = sorted(set(re.findall(r"VER:([\w.-]+)", raw.decode("utf-8-sig"))))
    return frame, pd.DataFrame(telemetry), pd.DataFrame(events), hashlib.sha256(raw).hexdigest()


def calibrate(frame: pd.DataFrame, fill_factor: float = FILL_FACTOR) -> pd.DataFrame:
    if not 0 < fill_factor <= 1:
        raise ValueError("Fill factor must be in (0, 1]")
    result = frame.copy()
    for raw_col, signed_col in (("c_raw", "c_signed"), ("v_raw", "v_signed"), ("p_raw", "p_signed")):
        result[signed_col] = result[raw_col].map(signed24)
    result["isc_A"] = result.c_signed / ADC_MAX * VREF_V / CURRENT_GAIN / SHUNT_OHM
    result["voc_V"] = result.v_signed / ADC_MAX * VREF_V * 2
    result["irradiance_W_m2"] = result.p_signed / ADC_MAX * VREF_V * 1000 / PYRANOMETER_MV_PER_W_M2
    result["assumed_fill_factor"] = fill_factor
    result["adc_clipped"] = result[["c_signed", "v_signed", "p_signed"]].abs().ge(ADC_MAX).any(axis=1)
    result["near_current_rail"] = result.c_signed.abs().ge(0.99 * ADC_MAX)
    result["invalid_polarity"] = result[["c_signed", "v_signed", "p_signed"]].lt(0).any(axis=1)
    result["illuminated"] = result.irradiance_W_m2 > MIN_IRRADIANCE_W_M2
    result["accepted"] = ~(result.adc_clipped | result.near_current_rail | result.invalid_polarity)
    result["quality_flag"] = np.select(
        [result.adc_clipped, result.near_current_rail, result.invalid_polarity],
        ["adc_clipped", "near_current_rail", "invalid_polarity"], default="accepted")
    result["estimated_mpp_W"] = (result.isc_A * result.voc_V * fill_factor).where(result.accepted)
    valid_response = result.accepted & result.illuminated
    result["instant_aperture_efficiency_pct"] = (
        100 * result.estimated_mpp_W / (APERTURE_M2 * result.irradiance_W_m2)
    ).where(valid_response)
    result["responsivity_A_m2_W"] = (result.isc_A / result.irradiance_W_m2).where(valid_response)
    return result


def accepted_edges(frame: pd.DataFrame, selected=None, max_gap_s: float = MAX_GAP_S):
    """Return a boolean mask for adjacent accepted sample pairs only."""
    n = len(frame)
    keep = np.ones(n, dtype=bool) if selected is None else np.asarray(selected, dtype=bool)
    dt = frame.timestamp.diff().dt.total_seconds().to_numpy()
    reboot = frame.uptime_s.diff().lt(0).to_numpy()
    edge = np.zeros(n, dtype=bool)
    edge[1:] = (keep[1:] & keep[:-1] & frame.accepted.to_numpy()[1:] & frame.accepted.to_numpy()[:-1]
                & frame.illuminated.to_numpy()[1:] & frame.illuminated.to_numpy()[:-1]
                & (dt[1:] > 0) & (dt[1:] <= max_gap_s) & ~reboot[1:])
    return edge


def integrate_edges(frame: pd.DataFrame, values, edge_mask) -> float:
    values = np.asarray(values, dtype=float)
    seconds = frame.timestamp.diff().dt.total_seconds().to_numpy()
    right = np.flatnonzero(edge_mask)
    left = right - 1
    return float(np.sum((values[left] + values[right]) * 0.5 * seconds[right]) / 3600)


def metrics(frame: pd.DataFrame, start=None, end=None, max_gap_s=MAX_GAP_S):
    t = frame.timestamp
    selected = np.ones(len(frame), dtype=bool)
    if start is not None:
        selected &= (t >= pd.Timestamp(start)).to_numpy()
    if end is not None:
        selected &= (t < pd.Timestamp(end)).to_numpy()
    edges = accepted_edges(frame, selected, max_gap_s)
    accepted = selected & frame.accepted.to_numpy() & frame.illuminated.to_numpy()
    obs = frame.loc[accepted]
    support_s = float(frame.timestamp.diff().dt.total_seconds().to_numpy()[edges].sum())
    elapsed_s = ((min(t.max(), pd.Timestamp(end)) if end is not None else t.max()) -
                 (max(t.min(), pd.Timestamp(start)) if start is not None else t.min())).total_seconds()
    energy = integrate_edges(frame, frame.estimated_mpp_W, edges)
    irradiation = integrate_edges(frame, frame.irradiance_W_m2, edges)
    ref_resp = float(frame.loc[frame.accepted & frame.illuminated, "responsivity_A_m2_W"].max())
    p99_resp = float(frame.loc[frame.accepted & frame.illuminated, "responsivity_A_m2_W"].quantile(.99))
    # Reference powers use observed Voc and the same common irradiance support.
    max_ref_power = ref_resp * frame.irradiance_W_m2 * frame.voc_V * frame.assumed_fill_factor
    p99_ref_power = p99_resp * frame.irradiance_W_m2 * frame.voc_V * frame.assumed_fill_factor
    max_ref_energy = integrate_edges(frame, max_ref_power, edges)
    p99_ref_energy = integrate_edges(frame, p99_ref_power, edges)
    peak = obs.loc[obs.estimated_mpp_W.idxmax()] if not obs.empty else None
    peak_isc = obs.loc[obs.isc_A.idxmax()] if not obs.empty else None
    peak_g = obs.loc[obs.irradiance_W_m2.idxmax()] if not obs.empty else None
    si_19 = irradiation * APERTURE_M2 * .19
    si_21 = irradiation * APERTURE_M2 * .21
    return {
        "start": str(obs.timestamp.min()) if not obs.empty else "",
        "end": str(obs.timestamp.max()) if not obs.empty else "",
        "sample_count": int(selected.sum()), "accepted_illuminated_samples": int(accepted.sum()),
        "integrated_intervals": int(edges.sum()), "elapsed_span_h": elapsed_s / 3600,
        "integrated_support_h": support_s / 3600,
        "coverage_pct": 100 * support_s / elapsed_s if elapsed_s > 0 else float("nan"),
        "irradiation_Wh_m2": irradiation, "incident_aperture_energy_Wh": irradiation * APERTURE_M2,
        "estimated_mpp_energy_Wh": energy,
        "energy_weighted_aperture_efficiency_pct": 100 * energy / (APERTURE_M2 * irradiation) if irradiation else float("nan"),
        "median_instant_aperture_efficiency_pct": float(obs.instant_aperture_efficiency_pct.median()) if not obs.empty else float("nan"),
        "peak_instant_aperture_efficiency_pct": float(obs.instant_aperture_efficiency_pct.max()) if not obs.empty else float("nan"),
        "empirical_max_reference_A_m2_W": ref_resp, "empirical_p99_reference_A_m2_W": p99_resp,
        "max_reference_ratio_pct": 100 * energy / max_ref_energy if max_ref_energy else float("nan"),
        "p99_reference_ratio_pct": 100 * energy / p99_ref_energy if p99_ref_energy else float("nan"),
        "silicon_19_energy_Wh": si_19, "silicon_21_energy_Wh": si_21,
        "cpv_to_silicon_19_ratio_pct": 100 * energy / si_19 if si_19 else float("nan"),
        "cpv_to_silicon_21_ratio_pct": 100 * energy / si_21 if si_21 else float("nan"),
        "peak_power_timestamp": str(peak.timestamp) if peak is not None else "",
        "peak_power_W": float(peak.estimated_mpp_W) if peak is not None else float("nan"),
        "aperture_efficiency_at_peak_power_pct": float(peak.instant_aperture_efficiency_pct) if peak is not None else float("nan"),
        "peak_power_isc_A": float(peak.isc_A) if peak is not None else float("nan"),
        "peak_power_voc_V": float(peak.voc_V) if peak is not None else float("nan"),
        "peak_power_irradiance_W_m2": float(peak.irradiance_W_m2) if peak is not None else float("nan"),
        "peak_isc_timestamp": str(peak_isc.timestamp) if peak_isc is not None else "",
        "peak_isc_A": float(peak_isc.isc_A) if peak_isc is not None else float("nan"),
        "peak_irradiance_timestamp": str(peak_g.timestamp) if peak_g is not None else "",
        "peak_irradiance_W_m2": float(peak_g.irradiance_W_m2) if peak_g is not None else float("nan"),
    }


def suncalc_frame():
    day = pd.Timestamp("2026-11-04")
    return pd.DataFrame([{"timestamp": day + pd.Timedelta(hours=int(t[:2]), minutes=int(t[3:])),
                          "azimuth_deg": az, "solar_elevation_deg": el,
                          "latitude_deg": 21.87304, "longitude_deg": -102.28282,
                          "source": "SunCalc PDF, printed local time; 2026-11-04"}
                         for t, az, el in SUNCALC_ROWS])


def make_events(events: pd.DataFrame, frame: pd.DataFrame, telemetry: pd.DataFrame):
    relevant = events[events.message.str.contains(
        r"\[CLOUD\]|F1:|F5:|Standby|Auto-arranque|Ejecutando comando|VER:", regex=True, na=False)].copy()
    relevant["event_type"] = np.select(
        [relevant.message.str.contains("Set Coords", na=False),
         relevant.message.str.contains("[CLOUD]", regex=False, na=False),
         relevant.message.str.contains(r"F1:|F5:", regex=True, na=False)],
        ["manual_coordinate_command", "remote_parameter_or_status_command", "controller_event"], default="logged_event")
    # Add algorithmic gaps with exact endpoints and source-line references.
    dt = frame.timestamp.diff().dt.total_seconds()
    gaps = []
    for i in np.flatnonzero(dt.to_numpy() > MAX_GAP_S):
        gaps.append({"timestamp": frame.timestamp.iloc[i], "message": f"Sampling gap: {dt.iloc[i]:.0f} s; no bridge in integration",
                     "source_line": int(frame.source_line.iloc[i]), "event_type": "sampling_gap",
                     "previous_source_line": int(frame.source_line.iloc[i - 1]),
                     "uptime_decreased": bool(frame.uptime_s.iloc[i] < frame.uptime_s.iloc[i - 1])})
    if not telemetry.empty:
        tel_dt = telemetry.timestamp.diff().dt.total_seconds()
        for i in np.flatnonzero(tel_dt.to_numpy() > MAX_GAP_S):
            gaps.append({"timestamp": telemetry.timestamp.iloc[i], "message": f"TEL telemetry gap: {tel_dt.iloc[i]:.0f} s",
                         "source_line": int(telemetry.source_line.iloc[i]), "event_type": "tel_sampling_gap",
                         "previous_source_line": int(telemetry.source_line.iloc[i - 1]),
                         "uptime_decreased": bool(telemetry.uptime_s.iloc[i] < telemetry.uptime_s.iloc[i - 1])})
        for i in np.flatnonzero(telemetry.uptime_s.diff().lt(0).to_numpy()):
            gaps.append({"timestamp": telemetry.timestamp.iloc[i], "message": "TEL uptime decreased",
                         "source_line": int(telemetry.source_line.iloc[i]), "event_type": "tel_uptime_decrease",
                         "previous_source_line": int(telemetry.source_line.iloc[i - 1]),
                         "uptime_decreased": True})
    relevant = relevant.assign(previous_source_line=np.nan, uptime_decreased=np.nan)
    return pd.concat([relevant, pd.DataFrame(gaps)], ignore_index=True, sort=False).sort_values("timestamp")


def add_solar_context(frame: pd.DataFrame, sun: pd.DataFrame):
    # Linear interpolation is only for plotting context; no extrapolation outside table.
    result = frame.copy()
    x = sun.timestamp.astype("int64").to_numpy()
    tq = result.timestamp.astype("int64").to_numpy()
    for col in ("azimuth_deg", "solar_elevation_deg"):
        result[col] = np.interp(tq, x, sun[col].to_numpy(), left=np.nan, right=np.nan)
        result.loc[(tq < x[0]) | (tq > x[-1]), col] = np.nan
    return result


def gap_breaks(frame):
    dt = frame.timestamp.diff().dt.total_seconds().to_numpy()
    return np.r_[False, (dt[1:] > MAX_GAP_S)] | ~frame.accepted.to_numpy()


def plot_outputs(frame, events, sun, summary):
    OUT.mkdir(parents=True, exist_ok=True)
    plt.rcParams.update({"font.family": "DejaVu Sans", "font.size": 9,
                         "axes.spines.top": False, "axes.spines.right": False,
                         "pdf.fonttype": 42, "ps.fonttype": 42})
    time = frame.timestamp
    breaks = gap_breaks(frame)
    fig, ax = plt.subplots(5, 1, figsize=(10, 10), sharex=True, layout="constrained")
    ax[0].plot(time, frame.irradiance_W_m2.mask(breaks), color="#d09128", lw=.65)
    ax[0].set_ylabel("Aperture-plane\nirradiance [W m$^{-2}$]")
    ax[1].plot(time, (1000 * frame.isc_A).mask(breaks), color="#167d8d", lw=.65)
    ax[1].set_ylabel("Short-circuit\ncurrent [mA]")
    ax[2].plot(time, frame.voc_V.mask(breaks), color="#7356a5", lw=.65)
    ax[2].set_ylabel("Open-circuit\nvoltage [V]")
    ax[3].plot(time, frame.estimated_mpp_W.mask(breaks), color="#167d8d", lw=.65, label="Estimated MPP power")
    ax[3].set_ylabel("Estimated power\n[W]")
    ax2 = ax[3].twinx(); ax2.plot(time, frame.instant_aperture_efficiency_pct.mask(breaks), color="#27824e", lw=.55, alpha=.8)
    ax2.set_ylabel("Aperture estimate [%]", color="#27824e")
    ax[4].plot(time, frame.x_steps.mask(breaks), color="#285b9a", lw=.6, label="X (azimuth axis)")
    ax[4].plot(time, frame.y_steps.mask(breaks), color="#c05b43", lw=.6, label="Y (elevation axis)")
    ax[4].set_ylabel("Logged coordinate\n[steps]"); ax[4].legend(ncol=2, frameon=False)
    for a in ax:
        a.grid(alpha=.2); a.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
        for g in pd.to_datetime(summary["gap_starts"]): a.axvline(g, color="#555555", lw=.6, ls=":", alpha=.6)
    for ts in pd.to_datetime(summary["command_times"]):
        for a in ax: a.axvline(ts, color="#ad3e38", lw=.55, ls="--", alpha=.5)
    ax[-1].set_xlabel("Local log clock (time-zone offset not encoded)")
    fig.suptitle("4 November 2026 field record: complete acquired span", x=.02, ha="left", fontsize=12)
    save_figure(fig, "field_day_overview")

    fig, ax = plt.subplots(2, 1, figsize=(9, 6), sharex=True, layout="constrained")
    ax[0].plot(sun.timestamp, sun.azimuth_deg, "o-", color="#3468a8", ms=2.5, lw=1)
    ax[0].set_ylabel("Sun azimuth [deg]")
    ax[1].plot(sun.timestamp, sun.solar_elevation_deg, "o-", color="#d27c25", ms=2.5, lw=1)
    ax[1].set_ylabel("Solar elevation [deg]"); ax[1].set_xlabel("SunCalc printed local time (15-min entries)")
    for a in ax:
        a.axvspan(frame.timestamp.min(), frame.timestamp.max(), color="#555555", alpha=.08, label="Acquired interval")
        a.grid(alpha=.2); a.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    ax[0].set_title("SunCalc context, Aguascalientes (21.87304° N, 102.28282° W)", loc="left")
    fig.text(.01, -.01, "Solar direction is not measured receiver pointing error or collector incidence.", fontsize=8)
    save_figure(fig, "solar_position_context")

    ctx = add_solar_context(frame, sun)
    colors = plt.colormaps["viridis"].with_extremes(over="#ad3e38")
    color_limit = summary["reference_p99"] * 1e3
    fig, ax = plt.subplots(2, 1, figsize=(10, 7), sharex=True, layout="constrained")
    sc = ax[0].scatter(ctx.timestamp, ctx.azimuth_deg, c=ctx.responsivity_A_m2_W * 1e3,
                       cmap=colors, vmin=0, vmax=color_limit, s=7, rasterized=True)
    ax[0].set_ylabel("Sun azimuth [deg]")
    cb = fig.colorbar(sc, ax=ax, pad=.01, extend="max"); cb.set_label("$I_{sc}/G$ [mA m$^2$ W$^{-1}$]")
    ax[1].scatter(ctx.timestamp, ctx.solar_elevation_deg, c=ctx.responsivity_A_m2_W * 1e3,
                   cmap=colors, vmin=0, vmax=color_limit, s=7, rasterized=True)
    ax[1].set_ylabel("Solar elevation [deg]"); ax[1].set_xlabel("Local log clock; SunCalc values linearly interpolated")
    for a in ax:
        a.grid(alpha=.2); a.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    ax[0].set_title("Solar position and measured normalized response", loc="left")
    fig.text(.01, -.01, "Ephemeris context only; not angular error. Red values exceed the P99 color scale; observations are retained.", fontsize=8)
    save_figure(fig, "solar_response_alignment")

    fig, ax = plt.subplots(2, 1, figsize=(9, 6), sharex=True, layout="constrained")
    response = frame.responsivity_A_m2_W * 1e3
    ax[0].plot(time, response.mask(breaks), color="#285b9a", lw=.65, label="$I_{sc}/G$")
    max_ref = summary["reference_max"] * 1e3; p99_ref = summary["reference_p99"] * 1e3
    for a in ax:
        a.axhline(max_ref, color="#b53d35", ls="--", lw=1, label="Full-record maximum reference")
        a.axhline(p99_ref, color="#27824e", ls="--", lw=1, label="Full-record P99 sensitivity")
    ax[1].plot(time, response.mask(breaks), color="#285b9a", lw=.65)
    anomaly = summary["anomaly_time"]
    ax[0].scatter([pd.Timestamp(anomaly)], [summary["anomaly_resp"] * 1e3], color="#b53d35", zorder=4, s=28)
    ax[0].set_ylim(0, max_ref * 1.08); ax[1].set_ylim(0, p99_ref * 1.35)
    ax[0].set_ylabel("Full range\n[mA m$^2$ W$^{-1}$]")
    ax[1].set_ylabel("Typical range\n[mA m$^2$ W$^{-1}$]"); ax[1].set_xlabel("Local log clock")
    ax[0].set_title("Irradiance-normalized response; full record retained", loc="left")
    for a in ax:
        for g in pd.to_datetime(summary["gap_starts"]): a.axvline(g, color="#555", lw=.6, ls=":")
        a.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M")); a.grid(alpha=.2)
    ax[0].legend(frameon=False, ncol=2, fontsize=8)
    save_figure(fig, "response_reference_sensitivity")

    fig, ax = plt.subplots(figsize=(8, 5), layout="constrained")
    labels = ["WiT-CPV estimated\nFF-based output", "Silicon scenario\n19% constant", "Silicon scenario\n21% constant"]
    vals = [summary["full"]["estimated_mpp_energy_Wh"], summary["full"]["silicon_19_energy_Wh"], summary["full"]["silicon_21_energy_Wh"]]
    colors = ["#167d8d", "#8798aa", "#b0bac5"]
    bars = ax.bar(labels, vals, color=colors, width=.62)
    ax.bar_label(bars, fmt="%.2f Wh", padding=3)
    ax.set_ylabel("Integrated scenario energy [Wh]"); ax.set_ylim(0, max(vals) * 1.2)
    ax.set_title("Equal-area scenarios in the measured aperture plane", loc="left")
    ax.text(.02, .96, "A = 0.002497 m²; coplanar irradiance; constant ηSi.\nNo temperature correction; not a measured panel comparison.",
            transform=ax.transAxes, va="top", fontsize=8)
    ax.grid(axis="y", alpha=.2)
    save_figure(fig, "fixed_silicon_scenario")


def save_figure(fig, stem):
    fig.savefig(OUT / f"{stem}.pdf", bbox_inches="tight")
    fig.savefig(OUT / f"{stem}.png", dpi=220, bbox_inches="tight")
    plt.close(fig)


def write_result_tex(summary, gap_text):
    m = summary["full"]
    macros = {
        "NovSamples": f"{summary['sample_count']:,}".replace(",", "{,}"),
        "NovHours": f"{m['integrated_support_h']:.3f}",
        "NovEnergy": f"{m['estimated_mpp_energy_Wh']:.3f}",
        "NovIrradiation": f"{m['irradiation_Wh_m2'] / 1000:.3f}",
        "NovEfficiency": f"{m['energy_weighted_aperture_efficiency_pct']:.2f}",
        "NovPeakPower": f"{m['peak_power_W']:.3f}",
        "NovApertureAtPeakPower": f"{m['aperture_efficiency_at_peak_power_pct']:.2f}",
        "NovPeakCurrent": f"{m['peak_power_isc_A'] * 1000:.1f}",
        "NovPeakVoltage": f"{m['peak_power_voc_V']:.3f}",
        "NovPeakIrradiance": f"{m['peak_power_irradiance_W_m2']:.1f}",
        "NovRelativeMax": f"{m['max_reference_ratio_pct']:.2f}",
        "NovRelativePercentile": f"{m['p99_reference_ratio_pct']:.2f}",
        "NovSiliconStandardRatio": f"{m['cpv_to_silicon_19_ratio_pct'] / 100:.3f}",
        "NovSiliconPremiumRatio": f"{m['cpv_to_silicon_21_ratio_pct'] / 100:.3f}",
        "NovSteadyEfficiency": f"{summary['phases'][2]['energy_weighted_aperture_efficiency_pct']:.2f}",
    }
    METRICS_TEX.write_text(
        "% Generated by analyze_2026_11_04.py.\n" + "\n".join(
            "\\newcommand{\\" + name + "}{" + value + "}" for name, value in macros.items()
        ) + "\n", encoding="utf-8")
    tex = rf"""% Generated by analyze_2026_11_04.py; do not hand-edit numerical values.
\subsection{{Acquisition coverage and operating chronology}}

The acquired span is 09:52:50--16:05:13, comprising \NovSamples{{}}
electrical/coordinate observations and {summary['tel_records']:,} separate
telemetry messages. All main observations pass the ADC quality checks and
exceed the analysis illumination threshold. There are no duplicate main
timestamps or decreasing main-sample uptimes. Five gaps exceed 15~s:
{gap_text}. Excluding these gaps gives \NovHours{{}}~h of integration
support out of {m['elapsed_span_h']:.3f}~h elapsed, or {m['coverage_pct']:.2f}\%
coverage. The audit retains every acquired observation; no clipped samples
or unobserved intervals are reconstructed.

The complete overview in Fig.~\ref{{fig:nov04_overview}} connects electrical
response with stage motion. There are {summary['cloud_commands']} logged remote
commands, including early tuning and repositioning. For example,
\texttt{{Delay=4}} was logged at 12:27:20, before the 12:32 boundary. Three event boundaries
divide the record without selecting a favorable performance threshold:
\texttt{{ErrMax=1100}} at 12:32:02, \texttt{{ErrMax=1000}} at 13:53:00,
and the first late manual-coordinate command at 15:12:42. All four sequential
phases are retained in Table~\ref{{tab:nov04_phases}}. No additional remote
commands are recorded inside the 12:32--13:53 phase; periodic controller
reference resets continue. This is an interval without logged operator
interventions, not a controlled comparison of firmware or settings.

\subsection{{Electrical output and aperture estimate}}

Applying Eqs.~\eqref{{eq:power_estimate}}--\eqref{{eq:aperture}} over the
complete accepted support yields \NovEnergy{{}}~Wh of estimated available-MPP
energy from \NovIrradiation{{}}~kWh\,m$^{{-2}}$ of measured aperture-plane irradiation,
with an energy-weighted aperture estimate of \NovEfficiency\%.
Table~\ref{{tab:nov04_hourly}} shows every acquired clock-hour bin, including
the initial low response and late decline, rather than replacing the complete
record with a midday-only integral.

At the maximum of the power estimate ({m['peak_power_timestamp']}), the
sequentially recorded values in that observation were
$I_{{sc}}={m['peak_power_isc_A']*1000:.1f}$ mA,
$V_{{oc}}={m['peak_power_voc_V']:.3f}$ V and $G={m['peak_power_irradiance_W_m2']:.1f}$ W m$^{{-2}}$,
giving {m['peak_power_W']:.3f}~W estimated output and an instantaneous aperture
estimate of \NovApertureAtPeakPower\% at that observation. Separate maxima are not combined
into a synthetic operating point. The 12:32--13:53 phase contributes
{summary['phases'][2]['estimated_mpp_energy_Wh']:.3f}~Wh over
{summary['phases'][2]['integrated_support_h']:.3f}~h, with an aperture estimate
of \NovSteadyEfficiency\%. The 15:00--16:00 clock-hour estimate falls to
{summary['hourly'][6]['energy_weighted_aperture_efficiency_pct']:.2f}\%, and
the phase after the first 15:12 manual command gives
{summary['phases'][4]['energy_weighted_aperture_efficiency_pct']:.2f}\%.
The record shows deterioration before and after that command, not evidence
that the command alone caused it.

\begin{{table}}[htbp]
\centering\small
\caption{{Complete acquisition and fixed clock-hour summaries. $H$ is aperture-plane
irradiation, $\widehat E$ is FF-based available energy and
$\widehat\eta_{{\mathrm{{ap}},E}}$ is the energy-weighted aperture estimate.
Bins contain only acquired observations. Both edge endpoints must lie inside a
bin, so subdivided integrals omit boundary-crossing edges and need not sum
exactly to the full-record integral.}}
\label{{tab:nov04_hourly}}
\begin{{tabular}}{{lrrrrr}}
\toprule
Window & $N$ & Support [h] & $H$ [kWh\,m$^{{-2}}$] & $\widehat E$ [Wh] & $\widehat\eta_{{\mathrm{{ap}},E}}$ [\%] \\
\midrule
{summary['table_rows']}
\bottomrule
\end{{tabular}}
\end{{table}}

\begin{{table}}[htbp]
\centering\small
\caption{{All sequential event-defined phases. Internal command boundaries
belong to the following phase; the final acquired sample is retained.
As in Table~\ref{{tab:nov04_hourly}},
edges crossing a boundary are omitted. Phases describe operation; they are not
matched experimental treatments.}}
\label{{tab:nov04_phases}}
\begin{{tabular}}{{lrrrr}}
\toprule
Phase boundary clocks & $N$ & Support [h] & $\widehat E$ [Wh] & $\widehat\eta_{{\mathrm{{ap}},E}}$ [\%] \\
\midrule
{summary['phase_table_rows']}
\bottomrule
\end{{tabular}}
\end{{table}}

\subsection{{Relative response and solar-position context}}

An empirical same-record $I_{{sc}}/G$ reference is highly sensitive to one point:
the maximum is {summary['reference_max']:.6f} A m$^2$ W$^{{-1}}$ at
{summary['anomaly_time']}, where $G={summary['anomaly_G']:.1f}$ W m$^{{-2}}$ and
$I_{{sc}}={summary['anomaly_I']*1000:.1f}$ mA; the P99 reference is
{summary['reference_p99']:.6f} A m$^2$ W$^{{-1}}$. The anomalous sample is retained.
In its four-minute neighborhood, median irradiance is
{summary['anomaly_neighborhood']['irradiance_median']:.1f} W m$^{{-2}}$
({summary['anomaly_neighborhood']['count']} samples), consistent with an isolated
irradiance dip rather than sustained low irradiance; the cause remains unresolved.
Across the accepted illuminated samples, median $I_{{sc}}/G$ is
{summary['responsivity_distribution']['median']:.6f} A m$^2$ W$^{{-1}}$
(interquartile range {summary['responsivity_distribution']['q25']:.6f}--
{summary['responsivity_distribution']['q75']:.6f}).
Full-record integrated relative-yield indicators are
{m['max_reference_ratio_pct']:.1f}\% (maximum reference) and
{m['p99_reference_ratio_pct']:.1f}\% (P99 sensitivity). These are internal,
reference-dependent indicators rather than independently validated tracking
efficiencies. Figure~\ref{{fig:nov04_response}} displays both the full
responsivity range and a typical-range view; the anomalous observation remains
in the analysis.

Interpolated SunCalc positions at the first and last observations are
{summary['solar_context']['first_azimuth_deg']:.2f}$^\circ$ and
{summary['solar_context']['last_azimuth_deg']:.2f}$^\circ$ azimuth, with elevations
{summary['solar_context']['first_elevation_deg']:.2f}$^\circ$ and
{summary['solar_context']['last_elevation_deg']:.2f}$^\circ$, respectively.
The printed table reaches {summary['solar_context']['maximum_elevation_deg']:.2f}$^\circ$
elevation at 12:30. Figure~\ref{{fig:nov04_solar_response}} relates this
trajectory to the measured normalized response. Logged coordinates span
$X={summary['coordinates']['x_min']:,}$--${summary['coordinates']['x_max']:,}$
and $Y={summary['coordinates']['y_min']:,}$--${summary['coordinates']['y_max']:,}$
steps. These are observed solar and coordinate ranges, not independently
validated acceptance angles or rated mechanical limits.

\subsection{{Conditional fixed-silicon reference}}
The coplanar, equal-area scenarios of Eq.~\eqref{{eq:silicon_scenario}}
yield
{m['silicon_19_energy_Wh']:.3f} and {m['silicon_21_energy_Wh']:.3f} Wh, respectively;
the corresponding estimated-CPV-to-scenario ratios are
\NovSiliconStandardRatio{{}} and \NovSiliconPremiumRatio{{}}.
Figure~\ref{{fig:nov04_silicon}} illustrates these assumptions, not measured
comparative yields. Changing the assumed CPV fill factor from 0.82 to 0.70 or
0.90 would change the full-record energy estimate to
{m['estimated_mpp_energy_Wh']*.70/FILL_FACTOR:.3f} or
{m['estimated_mpp_energy_Wh']*.90/FILL_FACTOR:.3f}~Wh, respectively.
This sensitivity range is illustrative, not a measured uncertainty interval.

\begin{{figure}}[htbp]\centering
\includegraphics[width=0.98\linewidth]{{../figures/2026-11-04/field_day_overview.pdf}}
\caption{{Complete 4 November acquisition. Vertical dotted lines mark sampling
gaps; dashed lines mark the three selected phase-boundary commands. Raw trends are not smoothed, and gaps are
not bridged. Electrical output is estimated using $FF=0.82$.}}
\label{{fig:nov04_overview}}\end{{figure}}

\begin{{figure}}[htbp]\centering
\includegraphics[width=0.92\linewidth]{{../figures/2026-11-04/solar_response_alignment.pdf}}
\caption{{SunCalc azimuth/elevation context aligned to the measured response.
SunCalc entries are spaced 15 minutes; intermediate positions are linearly
interpolated for plotting. Color encodes $I_{{sc}}/G$, not pointing error;
values above P99 remain plotted in the red over-range color.}}
\label{{fig:nov04_solar_response}}\end{{figure}}

\begin{{figure}}[htbp]\centering
\includegraphics[width=0.88\linewidth]{{../figures/2026-11-04/response_reference_sensitivity.pdf}}
\caption{{Irradiance-normalized response for the full record. The maximum
reference is dominated by a retained anomalous observation; P99 is shown only as
a sensitivity analysis. Neither line is an independent ideal-tracking standard.}}
\label{{fig:nov04_response}}\end{{figure}}

\begin{{figure}}[htbp]\centering
\includegraphics[width=0.75\linewidth]{{../figures/2026-11-04/fixed_silicon_scenario.pdf}}
\caption{{Conditional equal-area scenarios using the measured aperture-plane irradiance.
Silicon efficiency is held constant at 19\% or 21\%; no panel,
temperature correction or measured head-to-head comparison is represented.}}
\label{{fig:nov04_silicon}}\end{{figure}}

"""
    RESULT_TEX.write_text(tex, encoding="utf-8")


def run():
    frame, telemetry, events, source_hash = parse_log()
    frame = calibrate(frame)
    sun = suncalc_frame()
    frame = add_solar_context(frame, sun)
    frame["integrate_from_previous"] = accepted_edges(frame)
    frame["interval_from_previous_s"] = frame.timestamp.diff().dt.total_seconds()
    event_table = make_events(events, frame, telemetry)
    OUT.mkdir(parents=True, exist_ok=True)
    summary_full = metrics(frame)
    gaps = event_table[event_table.event_type == "sampling_gap"]
    commands = events[events.message.str.contains(r"\[CLOUD\]", regex=True, na=False)]
    # Phases use logged cut points, not a performance-selected window.
    milestones = []
    for message in ("[CLOUD] ErrMax=1100", "[CLOUD] ErrMax=1000",
                    "[CLOUD] Set Coords X:115000, Y:79000"):
        matches = events.loc[events.message.eq(message)]
        if len(matches) != 1:
            raise ValueError(f"Expected exactly one phase-boundary event: {message}")
        milestones.append(matches.iloc[0])
    boundaries = [event.timestamp for event in milestones]
    end = frame.timestamp.max() + pd.Timedelta(seconds=1)
    phase_defs = [("Full acquired record", frame.timestamp.min(), end),
                  ("Initial tuning and repositioning", frame.timestamp.min(), boundaries[0]),
                  ("After 12:32 ErrMax update", boundaries[0], boundaries[1]),
                  ("After 13:53 ErrMax update", boundaries[1], boundaries[2]),
                  ("After 15:12 manual coordinates", boundaries[2], end)]
    rows = []
    for name, start, end in phase_defs:
        m = metrics(frame, start, end)
        m["window"] = name
        m["window_start"] = str(start)
        m["window_end"] = str(min(end, frame.timestamp.max()))
        rows.append(m)
    # Fixed clock-hour bins retain late-session deterioration and expose support.
    hourly = []
    for h in range(9, 17):
        start = pd.Timestamp(f"2026-11-04 {h:02d}:00:00")
        m = metrics(frame, start, start + pd.Timedelta(hours=1))
        m["window"] = f"{h:02d}:00--{h+1:02d}:00"
        hourly.append(m)
    # Common thresholds and gap policies are explicit sensitivity checks.
    sensitivity = []
    for gap_limit in (10, 15, 30, 60):
        for light_limit in (0, 50, 100):
            saved = frame["illuminated"].copy()
            frame["illuminated"] = frame.irradiance_W_m2 > light_limit
            m = metrics(frame, max_gap_s=gap_limit)
            sensitivity.append({"max_gap_s": gap_limit, "minimum_irradiance_W_m2": light_limit,
                                "support_h": m["integrated_support_h"], "energy_Wh": m["estimated_mpp_energy_Wh"],
                                "aperture_efficiency_pct": m["energy_weighted_aperture_efficiency_pct"]})
            frame["illuminated"] = saved
    reference = frame.loc[frame.accepted & frame.illuminated, "responsivity_A_m2_W"]
    anomaly_idx = frame.responsivity_A_m2_W.idxmax()
    anomaly = frame.loc[anomaly_idx]
    anomaly_neighbor = frame.loc[(frame.timestamp >= anomaly.timestamp - pd.Timedelta(minutes=2)) &
                                 (frame.timestamp <= anomaly.timestamp + pd.Timedelta(minutes=2))]
    anomaly_neighbor_stats = {"count": len(anomaly_neighbor), "irradiance_min": float(anomaly_neighbor.irradiance_W_m2.min()),
                              "irradiance_median": float(anomaly_neighbor.irradiance_W_m2.median()),
                              "irradiance_max": float(anomaly_neighbor.irradiance_W_m2.max()),
                              "responsivity_median": float(anomaly_neighbor.responsivity_A_m2_W.median())}
    summary = {
        "source": str(SOURCE.relative_to(BASE)), "source_sha256": source_hash,
        "suncalc_source_sha256": hashlib.sha256(SUNCALC_SOURCE.read_bytes()).hexdigest(),
        "analysis_source_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        "source_bytes": SOURCE.stat().st_size, "source_lines": len(SOURCE.read_text(encoding="utf-8-sig").splitlines()),
        "sample_count": len(frame), "first_time": str(frame.timestamp.min()), "last_time": str(frame.timestamp.max()),
        "firmware": ", ".join(frame.attrs["firmware_versions"]), "cloud_commands": len(commands),
        "phase_boundary_events": [{"timestamp": str(event.timestamp), "message": event.message,
                                    "source_line": int(event.source_line)} for event in milestones],
        "tel_records": len(telemetry), "uptime_decreases": int(frame.uptime_s.diff().lt(0).sum()),
        "tel_uptime_decreases": int(telemetry.uptime_s.diff().lt(0).sum()),
        "tel_gaps_over_15s": int(telemetry.timestamp.diff().dt.total_seconds().gt(MAX_GAP_S).sum()),
        "tel_negative_step_deltas": int(telemetry.delta.lt(0).sum()),
        "duplicate_timestamps": int(frame.timestamp.duplicated().sum()),
        "gap_starts": [str(x) for x in gaps.timestamp],
        "gap_records": gaps[["timestamp", "message", "source_line", "previous_source_line", "uptime_decreased"]].to_dict("records"),
        "flags": {k: int(frame[k].sum()) for k in ("adc_clipped", "near_current_rail", "invalid_polarity", "illuminated", "accepted")},
        "full": summary_full, "phases": rows, "hourly": hourly,
        "solar_context": {"first_azimuth_deg": float(frame.azimuth_deg.iloc[0]),
                          "last_azimuth_deg": float(frame.azimuth_deg.iloc[-1]),
                          "first_elevation_deg": float(frame.solar_elevation_deg.iloc[0]),
                          "last_elevation_deg": float(frame.solar_elevation_deg.iloc[-1]),
                          "maximum_elevation_deg": float(sun.solar_elevation_deg.max()),
                          "clock_alignment": "Assumed log clock matches printed SunCalc local time; UTC offset unencoded"},
        "coordinates": {"x_min": int(frame.x_steps.min()), "x_max": int(frame.x_steps.max()),
                        "y_min": int(frame.y_steps.min()), "y_max": int(frame.y_steps.max())},
        "reference_max": float(reference.max()), "reference_p99": float(reference.quantile(.99)),
        "anomaly_time": str(anomaly.timestamp), "anomaly_line": int(anomaly.source_line),
        "anomaly_G": float(anomaly.irradiance_W_m2), "anomaly_I": float(anomaly.isc_A),
        "anomaly_resp": float(anomaly.responsivity_A_m2_W), "anomaly_neighborhood": anomaly_neighbor_stats,
        "responsivity_distribution": {
            "median": float(reference.median()), "q25": float(reference.quantile(.25)),
            "q75": float(reference.quantile(.75)), "p95": float(reference.quantile(.95)),
        },
        "calibration": {"adc_positive_full_scale": ADC_MAX, "vref_V": VREF_V, "current_gain": CURRENT_GAIN,
                        "shunt_ohm": SHUNT_OHM, "pyranometer_mV_per_W_m2": PYRANOMETER_MV_PER_W_M2,
                        "aperture_m2": APERTURE_M2, "fill_factor": FILL_FACTOR,
                         "max_gap_s": MAX_GAP_S, "minimum_irradiance_W_m2": MIN_IRRADIANCE_W_M2},
        "irradiance_geometry": {"pyranometer_coplanar_with_lens": True,
                                "provenance": "Experimenter confirmation on 2026-11-04",
                                "measured_quantity": "Global irradiance in the fixed lens-aperture plane",
                                "extra_cosine_factor_applied": False},
        "silicon_scenario": {"area_m2": APERTURE_M2, "efficiencies": [.19, .21],
                              "basis": "constant-efficiency equal-area coplanar silicon scenario with measured aperture-plane irradiance; no temperature correction; not measured panel"},
        "software": {"python": sys.version.split()[0], "numpy": np.__version__, "pandas": pd.__version__, "matplotlib": matplotlib.__version__},
    }
    # Text rows show complete record plus fixed clock bins without selecting a best interval.
    tab = [summary_full] + hourly
    summary["table_rows"] = "\n".join(
        f"{r['window'] if 'window' in r else 'Full acquired record'} & {r['accepted_illuminated_samples']} & "
        f"{r['integrated_support_h']:.3f} & {r['irradiation_Wh_m2'] / 1000:.3f} & {r['estimated_mpp_energy_Wh']:.4f} & "
        f"{r['energy_weighted_aperture_efficiency_pct']:.2f} " + "\\" * 2 for r in tab)
    phase_rows = rows[1:]
    summary["phase_table_rows"] = "\n".join(
        f"{pd.Timestamp(r['window_start']).strftime('%H:%M:%S')}--{pd.Timestamp(r['window_end']).strftime('%H:%M:%S')} & "
        f"{r['accepted_illuminated_samples']} & {r['integrated_support_h']:.3f} & "
        f"{r['estimated_mpp_energy_Wh']:.3f} & {r['energy_weighted_aperture_efficiency_pct']:.2f} " + "\\" * 2
        for r in phase_rows)
    # Mark command milestones relevant to phase boundaries without obscuring traces.
    summary["command_times"] = [str(t) for t in boundaries]
    summary["reference_max"] = float(reference.max()); summary["reference_p99"] = float(reference.quantile(.99))
    frame.to_csv(OUT / "audit_samples.csv", index=False, float_format="%.12g")
    event_table.to_csv(OUT / "event_chronology.csv", index=False)
    sun = suncalc_frame()
    sun.to_csv(OUT / "suncalc_2026-11-04.csv", index=False, float_format="%.5f")
    pd.DataFrame(rows + hourly).to_csv(OUT / "window_metrics.csv", index=False, float_format="%.10g")
    pd.DataFrame(sensitivity).to_csv(OUT / "sensitivity.csv", index=False, float_format="%.10g")
    (OUT / "analysis_summary.json").write_text(json.dumps(summary, indent=2, default=str, allow_nan=False) + "\n", encoding="utf-8")
    gap_text = "; ".join(
        f"{frame.loc[frame.source_line == int(g.previous_source_line), 'timestamp'].iloc[0].strftime('%H:%M:%S')}--"
        f"{pd.Timestamp(g.timestamp).strftime('%H:%M:%S')} ({int(float(g.message.split()[2]))} s)"
        for g in gaps.itertuples() if g.event_type == "sampling_gap")
    write_result_tex(summary, gap_text)
    plot_outputs(frame, event_table, sun, summary)
    return summary


if __name__ == "__main__":
    result = run()
    print(json.dumps({"samples": result["sample_count"], "full": result["full"],
                      "gaps": result["gap_records"], "source_sha256": result["source_sha256"]},
                     indent=2, default=str))