#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Estimacion de la irradiancia global REAL (correccion de la saturacion
del piranometro en 1024 W/m2) para una sesion de campo.

Ubicacion de las pruebas: Aguascalientes, Aguascalientes, Mexico.
    lat = +21.8756, lon = -102.2854, alt = 1896 m, TZ = UTC-6 (sin DST)

Metodo:
  1) Posicion solar por formulas estandar (declinacion de Cooper,
     ecuacion de tiempo de Spencer, correccion de longitud respecto al
     meridiano del huso UTC-6 a 90 W, masa de aire de Kasten-Young).
  2) Modelo de cielo despejado G_cs = A*exp(-B*AM) + C (transmisividad
     de Beer sobre la masa de aire + difuso de cielo claro), ajustado SOLO
     a muestras de cielo claro SIN saturar (el clipping a 1024 W/m2
     sesgaria el ajuste hacia abajo). Para anclar la forma de la curva,
     el ajuste se toma sobre el sobre-ensvolvente de la medida (muestra
     cercana al maximo movil de 30 min), porque las nubes solo bajan la
     curva de cielo despejado, nunca la suben.
  3) La correccion es MONOTONICA: a las muestras saturadas (G >= 1024)
     se les asigna G_est = max(1024, G_cs) (la correccion nunca baja la
     meda); al resto se les conserva la medida (nubes incluidas).
  4) Se recalculan con la misma ventana que el articulo (G > 50 W/m2, y
     G > 600 para las estadisticas de cielo claro): H, R_opt, P_ref,
     eta_tr, mediana P/Pref, fraccion P >= 0.95 P_ref y t95.

Uso:  py -3.12 clearsky.py <tag> [nofig]
Salida: articulo/resultados_<tag>/preview_irradiance_est.png
        articulo/resultados_<tag>/irradiance_est_<tag>.json
"""

import json
import sys
import tempfile
from pathlib import Path

import numpy as np
import pandas as pd

import matplotlib
matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt

import parse_log as pl

ROOT = Path(__file__).resolve().parent

# Aguascalientes, Ags.
LAT = 21.8756
LON_W = 102.2854     # longitud oeste positiva
ALT_M = 1896.0
L_STD_W = 90.0       # meridiano del huso UTC-6
G_CLIP = 1023.9      # ADC a tope => G = 1024.0 W/m2


def solar_geometry(ts: pd.Series):
    """Elevacion solar y masa de aire (Kasten-Young) a partir de la hora local."""
    n = (ts.dt.dayofyear.values * 86400.0
         + ts.dt.hour.values * 3600.0 + ts.dt.minute.values * 60.0
         + ts.dt.second.values) / 86400.0
    rad = np.pi / 180.0
    phi = LAT * rad
    B = 2 * np.pi / 365.0 * (n - 81.0)
    delta = 23.45 * rad * np.sin(2 * np.pi / 365.0 * (n + 284.0))
    eot_min = 229.18 * (0.000075 + 0.001868 * np.cos(B)
                         - 0.032077 * np.sin(B)
                         - 0.014615 * np.cos(2 * B)
                         - 0.040849 * np.sin(2 * B))
    clock_min = (ts.dt.hour.values * 3600 + ts.dt.minute.values * 60
                 + ts.dt.second.values) / 60.0
    l_corr = 4.0 * (LON_W - L_STD_W)          # minutos (+ si esta al oeste)
    hour_angle = 15.0 * ((clock_min - l_corr + eot_min) / 60.0 - 12.0)
    sin_el = (np.sin(phi) * np.sin(delta)
              + np.cos(phi) * np.cos(delta) * np.cos(hour_angle * rad))
    el = np.degrees(np.arcsin(np.clip(sin_el, -1, 1)))
    z = 90.0 - el
    m_air = np.where(z < 75.0,
                     1.0 / (np.cos(z * rad) + 0.14 * (93.85 - z) ** (-1.375)),
                     np.inf)
    return el, m_air


def fit_clearsky(el_deg, am, G, cloud_mask, t_idx):
    """Ajusta el modelo G_cs = A*exp(-B*AM) + C a las muestras claras no saturadas.

    Conjunto de ajuste: el sobre-ensvolvente de la medida (muestra >= 0.98
    del maximo movil de 30 min), fuera de transitorios nubosos, sin
    saturar (G < 1015) y con el sol por encima de 12 grados. Al anclarse
    en la parte alta de la curva, la extrapolacion al pico (que cae en la
    zona saturada del piranometro) es corta y confiable.
    """
    roll = pd.Series(G, index=pd.DatetimeIndex(t_idx)).rolling(
        "30min", center=True, min_periods=5).max().values
    fitm = ((el_deg > 12.0) & (G > 500.0) & (G < 1015.0)
            & (~cloud_mask) & (G >= 0.98 * roll))
    Gs = G[fitm]
    AMs = am[fitm]
    n_fit = int(fitm.sum())
    if n_fit < 20:
        raise RuntimeError("Muy pocas muestras claras para ajustar el modelo "
                           f"(n={n_fit}).")

    def sse_of(A, Bc, C):
        return float(np.sum((Gs - (A * np.exp(-Bc * AMs) + C)) ** 2))

    best = None
    for A in range(900, 1501, 50):
        for Bc in np.arange(0.05, 0.61, 0.02):
            for C in range(0, 301, 10):
                s = sse_of(A, Bc, C)
                if best is None or s < best[0]:
                    best = (s, float(A), float(Bc), float(C))
    # refinamiento local
    _, A0, B0, C0 = best
    for A in range(int(A0) - 50, int(A0) + 51, 10):
        for Bc in np.arange(B0 - 0.02, B0 + 0.021, 0.005):
            for C in range(int(C0) - 10, int(C0) + 11, 5):
                if A <= 0 or Bc <= 0 or C < 0:
                    continue
                s = sse_of(A, Bc, C)
                if best is not None and s < best[0]:
                    best = (s, float(A), float(Bc), float(C))
    sse, A, Bc, C = best
    rms = float(np.sqrt(sse / n_fit))
    am_finite = np.isfinite(am)
    am_safe = np.where(am_finite, am, 0.0)
    # La curva solo se evalua donde la masa de aire es finita (el>=15). Para
    # 0<el<15 la formula de Kasten-Young diverge (am=inf); reportar
    # A*exp(-B*0)=A ahi seria fisicamente incorrecto.
    curve = np.where((el_deg > 0.0) & am_finite,
                     A * np.exp(-Bc * am_safe) + C, 0.0)
    # Pico del modelo a mediodia (masa de aire minima) = irradiancia maxima
    # de cielo despejado en el sitio.
    if am_finite.any():
        am_min = float(np.min(am[am_finite]))
    else:
        am_min = 1.0
    model_peak = A * np.exp(-Bc * am_min) + C
    return A, Bc, C, rms, n_fit, curve, model_peak


def run(tag, nofig=False):
    from compare_sessions import normalize
    src = ROOT / f"{tag}.txt"
    if not src.exists():
        raise SystemExit(f"Log {src.name} no encontrado en la raiz.")
    tmp = Path(tempfile.gettempdir()) / (tag + "_clearsky.txt")
    normalize(src, tmp)
    df_s, df_e = pl.parse_log(tmp)
    df, ph = pl.build_physical(df_s, df_e)

    G = df["G_POA"].values
    Voc = df["Voc"].values
    Isc = df["Isc"].values
    Pc = df["P_CPV"].values
    t = df["timestamp"]

    el, am = solar_geometry(df["timestamp"])
    ces = pl.cloud_events(df)
    cloud_mask = np.zeros(len(df), dtype=bool)
    for c in ces:
        t0 = pd.Timestamp(c["t_start"])
        t1 = t0 + pd.Timedelta(seconds=c["dur_s"])
        cloud_mask |= (df["timestamp"] >= t0) & (df["timestamp"] <= t1)

    A, Bc, C, rms, n_fit, curve, model_peak = fit_clearsky(
        el, am, G, cloud_mask, t.values)
    clipped = G >= G_CLIP
    n_clip = int(clipped.sum())
    dt_med = float(np.median(df["dt_s"]))

    # correccion monotona: lo saturado toma el modelo (nunca por debajo
    # de los 1024 medidos); lo demas conserva la medida
    G_est = np.where(clipped, np.maximum(1024.0, curve), G)

    # R_opt y P_ref: con y sin correccion
    lit = G > 600
    R_now = float((Isc[lit] / G[lit]).max())
    P_ref_now = df["P_ref"].values
    litc = G_est > 600
    R_corr = float((Isc[litc] / G_est[litc]).max())
    P_ref_corr = R_corr * G_est * Voc * pl.FF

    # ventanas identicas a las del articulo (metrics.json)
    m_full = pl.window_metrics(df, "Full (G>50)")
    m_clear = pl.window_metrics(df[df["G_POA"] > 600], "Clear (G>600)")
    df_c = df.copy()
    df_c["G_POA"] = G_est
    df_c["P_ref"] = P_ref_corr
    m_full_c = pl.window_metrics(df_c, "Full (G>50)")
    m_clear_c = pl.window_metrics(df_c[df_c["G_POA"] > 600], "Clear (G>600)")

    df_mod = df.copy()
    df_mod["P_ref"] = P_ref_corr
    t95_now = pl.recovery_t95(df, df_e)
    t95_corr = pl.recovery_t95(df_mod, df_e)

    # MAE con la referencia corregida (ventana de cielo claro del articulo)
    Isc_ref_corr = R_corr * G_est
    wc = (G > 600)
    MAE_I_now = float(np.mean(np.abs((R_now * G) - Isc)) * 1000.0)
    MAE_I_corr = float(np.mean(np.abs(Isc_ref_corr - Isc)[wc])) * 1000.0
    MAE_P_now = m_clear["MAE_P_mW"]
    MAE_P_corr = float(np.mean(np.abs(P_ref_corr - Pc)[wc])) * 1000.0

    print("\n" + "=" * 72)
    print(f"   CORRECCION DE SATURACION (1024 W/m2)  -  {tag}")
    print(f"   Aguascalientes, Ags., Mexico (21.8756 N, 102.2854 W, 1896 m)")
    print("=" * 72)
    print(f"  Modelo cielo despejado: G_cs = {A:.0f} * exp(-{Bc:.3f} * AM) + {C:.0f}")
    print(f"    rms sobre n={n_fit} muestras claras no saturadas: {rms:.0f} W/m2")
    print(f"  Pico: medido = {G.max():.0f} W/m2 (tope del ADC) | "
          f"estimado = {G_est.max():.0f} W/m2")
    print(f"  Muestras saturadas: {n_clip} "
          f"({100.0 * n_clip / len(df):.1f} % del log, ~{n_clip * dt_med / 60:.0f} min)")
    print(f"  H (kWh/m2, G>50):    {m_full['H_kWh_m2']} -> "
          f"{m_full_c['H_kWh_m2']} "
          f"(+{100.0 * (m_full_c['H_kWh_m2'] / m_full['H_kWh_m2'] - 1):.1f} %)")
    print(f"  R_opt (uA/W/m2):     {R_now * 1e6:.0f} -> {R_corr * 1e6:.0f}")
    print(f"  eta_tr (G>50):       {m_full['eta_tr_pct']} % -> "
          f"{m_full_c['eta_tr_pct']} %")
    print(f"  eta_tr (G>600):      {m_clear['eta_tr_pct']} % -> "
          f"{m_clear_c['eta_tr_pct']} %")
    print(f"  MAE_I (G>600):       {MAE_I_now:.1f} mA -> "
          f"{MAE_I_corr:.1f} mA")
    print(f"  MAE_P (G>600):       {MAE_P_now} mW -> {MAE_P_corr:.1f} mW")
    print(f"  med P/Pref (G>600):  {m_clear['ratio_P_Pref_median'] if 'ratio_P_Pref_median' in m_clear else float(np.median(Pc[wc] / P_ref_now[wc])):.3f} -> "
          f"{float(np.median(Pc[wc] / P_ref_corr[wc])):.3f}")
    if len(t95_now) == len(t95_corr) and len(t95_now):
        print("\n  t95 (s):   antes -> corregido")
        for a, b in zip(t95_now.itertuples(), t95_corr.itertuples()):
            print(f"    {a.trigger:>14} @ {str(a.t_trigger)[11:19]}: "
                  f"{a.t95_s} -> {b.t95_s}")

    if not nofig:
        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 4.4),
                                        sharex=True)
        ax1.plot(t, Pc * 1000, color="tab:red", lw=0.8, label="P_CPV")
        ax1.plot(t, P_ref_now * 1000, color="k", lw=1.1,
                 label="P_ref (G medida)")
        ax1.plot(t, P_ref_corr * 1000, color="tab:blue", ls="--", lw=1.1,
                 label="P_ref (G corregida)")
        ax1.set_ylabel("Potencia (mW)")
        ax1.set_xlabel("Hora local (UTC-6)")
        ax1.legend(fontsize=8)
        ax1.grid(alpha=0.3)

        if n_clip:
            idx = np.where(clipped)[0]
            start = 0
            for j in range(1, len(idx) + 1):
                if j == len(idx) or idx[j] != idx[j - 1] + 1:
                    ax2.axvspan(t.iloc[idx[start]], t.iloc[idx[j - 1]],
                                 color="tab:orange", alpha=0.18)
                    start = j
        ax2.plot(t, G, color="tab:orange", lw=0.8, label="G medida")
        ax2.plot(t, G_est, color="tab:orange", ls="--", lw=1.3,
                 label="G estimada (corregida)")
        ax2.plot(t, curve, color="gray", ls=":", lw=0.8,
                 label="Modelo cielo despejado")
        ax2.axhline(1024, color="gray", ls="-.", lw=0.8)
        ax2.text(t.iloc[1], 1030, "saturacion 1024 W/m2",
                 fontsize=8, color="dimgray")
        ax2.set_ylabel("Irradiancia global (W/m2)")
        ax2.set_xlabel("Hora local (UTC-6)")
        ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
        ax2.legend(fontsize=8, loc="upper right")
        ax2.grid(alpha=0.3)
        fig.tight_layout()

        outdir = ROOT / "articulo" / f"resultados_{tag}"
        outdir.mkdir(parents=True, exist_ok=True)
        fig.savefig(outdir / "preview_irradiance_est.png", dpi=150)

    res = {
        "tag": tag,
        "model": "G_cs = A*exp(-B*AM) + C, Ajustado al sobre-ensvolvente "
                 "claro no saturado",
        "model_scale_A_Wm2": round(A, 1), "model_B": round(Bc, 4), "model_C": C,
        "model_rms_clear_Wm2": round(rms, 1), "n_fit_clear": n_fit,
        "model_peak_Wm2": round(model_peak, 0),
        "G_max_measured_Wm2": round(float(G.max()), 1),
        "G_max_estimated_Wm2": round(float(G_est.max()), 1),
        "n_clipped": n_clip,
        "clipped_pct_log": round(100 * n_clip / len(df), 2),
        "H_kWh_m2_full_measured": m_full["H_kWh_m2"],
        "H_kWh_m2_full_corrected": m_full_c["H_kWh_m2"],
        "R_opt_uA_per_Wm2_measured": round(R_now * 1e6, 1),
        "R_opt_uA_per_Wm2_corrected": round(R_corr * 1e6, 1),
        "eta_tr_pct_full_measured": m_full["eta_tr_pct"],
        "eta_tr_pct_full_corrected": m_full_c["eta_tr_pct"],
        "eta_tr_pct_clear_measured": m_clear["eta_tr_pct"],
        "eta_tr_pct_clear_corrected": m_clear_c["eta_tr_pct"],
        "MAE_I_mA_clear_measured": round(MAE_I_now, 2),
        "MAE_I_mA_clear_corrected": round(MAE_I_corr, 2),
        "MAE_P_mW_clear_measured": MAE_P_now,
        "MAE_P_mW_clear_corrected": round(MAE_P_corr, 2),
        "ratio_P_Pref_median_clear_measured":
            round(float(np.median(Pc[wc] / P_ref_now[wc])), 3),
        "ratio_P_Pref_median_clear_corrected":
            round(float(np.median(Pc[wc] / P_ref_corr[wc])), 3),
        "location": "Aguascalientes, Ags., Mexico "
                    "(21.8756 N, 102.2854 W, 1896 m, UTC-6)",
    }
    outdir = ROOT / "articulo" / f"resultados_{tag}"
    outdir.mkdir(parents=True, exist_ok=True)
    (outdir / f"irradiance_est_{tag}.json").write_text(
        json.dumps(res, indent=2, ensure_ascii=False), encoding="utf-8")
    print(f"\n[ok] JSON (y figura) en {outdir}")


if __name__ == "__main__":
    tag = sys.argv[1] if len(sys.argv) > 1 else "21092025"
    nofig = len(sys.argv) > 2 and sys.argv[2] == "nofig"
    run(tag, nofig=nofig)