simulator.py

"""
Physics-based PV simulator using the simplified single-diode model.

Single-diode model (simplified, no numerical solving):
    I = Iph - I0*(exp((V+I*Rs)/(n*Vt)) - 1) - (V+I*Rs)/Rsh

At the maximum power point, observable outputs are:
    Isc  — short-circuit current  (linear in G, slight +ve temp coefficient)
    Voc  — open-circuit voltage   (logarithmic in G, strong -ve temp coefficient)
    FF   — fill factor            (slight -ve temp coefficient via Rs heating)
    Pmax — Isc * Voc * FF

System: 100kWp c-Si, 15 strings x 22 modules in series (see config.py)
    V_OC_SYS = 22 x 40V = 880V    I_SC_SYS = 15 x 9.8A = 147A
    P_STC    = 880 x 147 x 0.78 = 100.9kWp

Fault physics:
    string_anomaly     — 3 strings at 15% Isc in afternoon (shading/cell mismatch)
                         Isc drops ~17%, Voc unchanged (parallel strings hold voltage)
    bypass_diode       — 8 modules in one string with one shorted bypass diode each
                         1/3 of each module's cells permanently bypassed
                         Voc of that string: 22*40 - 8*(40/3) = 773V vs 880V normal
                         System Voc: (14*880 + 773)/15 = 869V -> ~1.3% Voc drop
                         Isc slightly elevated (lower Rs in bypassed section)
                         FF degrades (current mismatch increases effective Rs)
    module_fault       — 3 strings open-circuit (combiner box or severe cell failure)
                         Isc drops to 12/15 = 80%, Voc slightly lower
    gradual_degradation — EVA yellowing (Isc down), metallisation corrosion (FF down),
                          microcracking (Rsh down -> FF down), ~1.5%/yr

Weather sources (priority):
    1. PVGIS ERA5 historical irradiance + ambient temp (past dates, interpolated)
    2. Monthly-calibrated synthetic noise from PVGIS SD data (today/future)
"""

from datetime import date, datetime, timedelta, timezone
from functools import lru_cache

import numpy as np

import config as _cfg
from config import (
    SYSTEM_PEAK_KW, SYSTEM_UTC_OFFSET, SYSTEM_LAT, SYSTEM_LON,
    N_STRINGS, N_SERIES, N_BYPASS, V_OC_MOD, I_SC_MOD, FF_STC,
    ALPHA_ISC, BETA_VOC, NOCT, DIODE_N, CELLS_PER_MOD, SYSTEM_LOSS,
)
from pvgis_client import build_lookup_table, get_expected, get_historical_conditions

# ── Derived system-level STC parameters ──────────────────────────────────────
G_REF     = 1000.0                     # W/m²  STC irradiance
T_REF     = 25.0                       # °C    STC temperature
V_OC_SYS  = N_SERIES  * V_OC_MOD      # 880V  system open-circuit voltage
I_SC_SYS  = N_STRINGS * I_SC_MOD      # 147A  system short-circuit current
N_CELLS_S = N_SERIES  * CELLS_PER_MOD # 1320  total cells in series (Vt scaling)

# ── Monthly irradiance variability — from PVGIS SD data, Ahmedabad 100kWp ────
_MONTHLY_CV = {
    1: 0.14, 2: 0.18, 3: 0.15, 4: 0.09,
    5: 0.12, 6: 0.34, 7: 0.57, 8: 0.78,
    9: 0.66, 10: 0.23, 11: 0.15, 12: 0.13,
}

_PVGIS_HISTORICAL_YEAR = 2019   # ERA5 year used as weather proxy for past dates


# ── Physics functions ─────────────────────────────────────────────────────────

def _cell_temp(G: float, tamb: float) -> float:
    """
    NOCT model: Tcell = Tamb + (NOCT - 20) * G / 800
    NOCT = 45 deg C for standard free-standing c-Si modules.
    Replaces the old irradiance_fraction * 20 approximation.
    """
    return tamb + (NOCT - 20.0) * G / 800.0


def _thermal_voltage(tcell: float) -> float:
    """Thermal voltage Vt = kT/q (volts). k = 8.617e-5 eV/K."""
    return 8.617e-5 * (tcell + 273.15)


def _pv_output(G: float, tamb: float,
               isc_frac: float = 1.0,
               voc_frac: float = 1.0,
               ff_frac:  float = 1.0) -> dict:
    """
    Compute PV system output from irradiance and ambient temperature.

    Equations:
        Isc(G, T) = Isc_ref * (G/G_ref) * (1 + alpha_Isc * dT)
        Voc(G, T) = Voc_ref * (1 + beta_Voc * dT)
                    + N_cells * n * Vt * ln(G/G_ref)   [irradiance correction]
        FF(T)     = FF_ref * (1 - 0.0015 * max(dT, 0)) [Rs heating effect]
        Pmax      = Isc * Voc * FF

    Fault multipliers (isc_frac, voc_frac, ff_frac) are applied independently
    so the classifier can observe WHICH parameter deviated.
    """
    if G < 10.0:
        return {
            "power_w": 0, "voltage": 0.0, "current": 0.0,
            "fill_factor": FF_STC, "irradiance": round(G, 1),
            "cell_temp": round(tamb, 1),
        }

    tcell   = _cell_temp(G, tamb)
    dt      = tcell - T_REF
    g_ratio = G / G_REF

    # Isc: linear in irradiance, slight positive temperature coefficient
    isc = I_SC_SYS * g_ratio * (1.0 + ALPHA_ISC * dt) * isc_frac

    # Voc: negative temp coefficient + logarithmic irradiance dependence
    voc_t   = V_OC_SYS * (1.0 + BETA_VOC * dt)
    voc_irr = N_CELLS_S * DIODE_N * _thermal_voltage(tcell) * np.log(max(g_ratio, 1e-6))
    voc     = max((voc_t + voc_irr) * voc_frac, 0.0)

    # Fill factor: slight degradation at high cell temp (thermal Rs increase)
    ff = float(np.clip(FF_STC * (1.0 - 0.0015 * max(dt, 0.0)) * ff_frac, 0.50, FF_STC))

    # DC array power at the maximum power point, then apply system loss
    # (inverter, wiring, soiling, mismatch) to get AC power at the meter --
    # the same quantity PVGIS reports with loss=14, so measured and the PVGIS
    # baseline are directly comparable. Voltage/current below stay DC (sensor
    # values at the array), so the fault classifier is unaffected.
    pmax = isc * voc * ff * (1.0 - SYSTEM_LOSS)

    return {
        "power_w":     round(max(pmax, 0.0)),
        "voltage":     round(voc, 1),        # DC array Voc — fault-sensitive
        "current":     round(isc, 2),         # DC array Isc — fault-sensitive
        "fill_factor": round(ff, 4),          # FF — fault-sensitive
        "irradiance":  round(G, 1),           # W/m² plane-of-array
        "cell_temp":   round(tcell, 1),       # °C cell temperature via NOCT
    }


def _fault_factors(dt: datetime, fault_mode: str) -> tuple[float, float, float]:
    """
    Return (isc_frac, voc_frac, ff_frac) encoding the physical fault signature.

    Each fault modifies specific electrical parameters:
        string_anomaly:  Isc drops (current mismatch), Voc unchanged
        bypass_diode:    Voc drops (fewer cells in series), Isc slightly up, FF down
        module_fault:    Isc drops severely (strings lost), Voc slightly down
        degradation:     All slowly decline (Isc via EVA, FF via Rs/Rsh)
    """
    if fault_mode == "none":
        return 1.0, 1.0, 1.0

    if fault_mode == "string_anomaly":
        # 3 strings at 15% Isc in afternoon — geometry-driven shading/mismatch
        if dt.hour >= 11:
            isc_frac = (N_STRINGS - 3 + 3 * 0.15) / N_STRINGS  # (12 + 0.45)/15 = 0.830
            return isc_frac, 1.0, 1.0
        return 1.0, 1.0, 1.0

    if fault_mode == "bypass_diode":
        # 8 modules with one shorted diode (1/3 cells bypassed), all in one string
        voc_string_faulty = N_SERIES * V_OC_MOD - 8 * (V_OC_MOD / N_BYPASS)
        voc_sys_frac = ((N_STRINGS - 1) * V_OC_SYS + voc_string_faulty) / (N_STRINGS * V_OC_SYS)
        # Bypassed cells lower string Rs -> slight Isc increase; mismatch degrades FF
        return 1.02, voc_sys_frac, 0.96

    if fault_mode == "module_fault":
        # 3 strings open-circuit (combiner box failure or catastrophic cell damage)
        isc_frac = (N_STRINGS - 3) / N_STRINGS  # 12/15 = 0.800
        return isc_frac, 0.97, 0.95

    if fault_mode == "gradual_degradation":
        # EVA yellowing (Isc), metallisation corrosion + microcracking (FF)
        return 0.985, 0.997, 0.988

    if fault_mode == "ground_fault":
        # Insulation / ground leakage: a shunt path to ground lowers effective
        # Rsh -> moderate current loss and FF degradation, Voc roughly stable.
        return 0.93, 0.99, 0.93

    return 1.0, 1.0, 1.0


# ── Weather functions ─────────────────────────────────────────────────────────

@lru_cache(maxsize=400)
def _synthetic_day_irr_factor(day_ordinal: int, month: int) -> float:
    """
    Day-level irradiance multiplier for synthetic (non-PVGIS) dates.
    Seeded by ordinal for reproducibility.
    """
    rng  = np.random.default_rng(seed=day_ordinal)
    cv   = _MONTHLY_CV.get(month, 0.20)
    # Center a typical day on the TMY mean (1.0), not 0.85: the baseline this is
    # compared against IS the TMY mean (which already averages in cloudy days),
    # so a healthy day should track it. Variability comes from the monthly CV.
    return float(np.clip(rng.normal(1.0, cv), 0.35, 1.15))


def _cloud_noise(rng: np.random.Generator, month: int) -> float:
    """Sub-hourly cloud flicker, scaled to monthly variability."""
    cv  = _MONTHLY_CV.get(month, 0.20)
    std = float(np.clip(cv * 0.3, 0.05, 0.25))
    return float(np.clip(rng.normal(1.0, std), 0.0, 1.3))


def _get_slot_conditions(target: date, pvgis_lookup: dict,
                         hour: int, minute: int,
                         day_rng: np.random.Generator,
                         synth_irr_factor: float,
                         lat: float, lon: float) -> tuple[float, float]:
    """
    Return (irradiance_wm2, tamb_c) for a 5-min slot.

    Past dates: linearly interpolate PVGIS ERA5 hourly irradiance and temp
    between current and next hour, then add sub-hourly noise.
    Today/future: scale TMY mean by synthetic day factor and cloud noise.
    """
    month = target.month
    exp   = get_expected(month, hour, pvgis_lookup)

    if target < date.today():
        irr0, t0 = get_historical_conditions(
            month, target.day, hour, _PVGIS_HISTORICAL_YEAR, lat, lon)
        irr1, t1 = get_historical_conditions(
            month, target.day, (hour + 1) % 24, _PVGIS_HISTORICAL_YEAR, lat, lon)

        if irr0 is not None and irr1 is not None:
            frac  = minute / 60.0
            G     = max(irr0 + frac * (irr1 - irr0), 0.0)
            G    *= float(np.clip(day_rng.normal(1.0, 0.04), 0.0, 1.2))  # sub-hourly flicker
            tamb  = (t0 + frac * (t1 - t0)) if t1 is not None else t0
            return G, tamb

    # Synthetic fallback
    G_mean = exp.get("mean_irradiance_wm2", 0.0)
    tamb   = exp.get("mean_temp_c", 25.0) + float(day_rng.normal(0.0, 1.5))
    G      = G_mean * synth_irr_factor * _cloud_noise(day_rng, month)
    return max(G, 0.0), tamb


# ── Core reading generator ────────────────────────────────────────────────────

def _readings_for_date(target: date, pvgis_lookup: dict,
                       fault_mode: str = "none",
                       lat: float = SYSTEM_LAT,
                       lon: float = SYSTEM_LON) -> list[dict]:
    """
    Generate all 288 x 5-min readings for a date using the physics model.
    Each reading contains power, DC voltage (Voc), DC current (Isc),
    irradiance W/m², cell temperature, fill factor, and cumulative energy.
    """
    month        = target.month
    day_rng      = np.random.default_rng(seed=target.toordinal() + 1000)
    synth_f      = _synthetic_day_irr_factor(target.toordinal(), month)
    date_str     = target.strftime("%Y%m%d")
    readings     = []
    cumulative_wh = 0.0

    for slot in range(288):
        total_min = slot * 5
        hour      = total_min // 60
        minute    = total_min % 60
        dt        = datetime(target.year, month, target.day, hour, minute)

        G, tamb = _get_slot_conditions(
            target, pvgis_lookup, hour, minute, day_rng, synth_f, lat, lon)

        isc_frac, voc_frac, ff_frac = _fault_factors(dt, fault_mode)
        out = _pv_output(G, tamb, isc_frac, voc_frac, ff_frac)

        cumulative_wh += out["power_w"] * (5 / 60)

        readings.append({
            "date":        date_str,
            "time":        f"{hour:02d}:{minute:02d}",
            "energy_wh":   round(cumulative_wh),
            "power_w":     out["power_w"] if out["power_w"] > 0 else None,
            "temperature": out["cell_temp"],    # cell temp via NOCT model
            "voltage":     out["voltage"],       # DC array Voc (fault-sensitive)
            "current":     out["current"],       # DC array Isc (fault-sensitive)
            "irradiance":  out["irradiance"],    # W/m² plane-of-array
            "fill_factor": out["fill_factor"],   # fault-sensitive
        })

    return readings


# ── Module-level PVGIS lookup cache ──────────────────────────────────────────

_pvgis_lookup: dict | None = None

def _get_lookup() -> dict:
    global _pvgis_lookup
    if _pvgis_lookup is None:
        _pvgis_lookup = build_lookup_table()
    return _pvgis_lookup


# ── Public API — same interface as before ─────────────────────────────────────

def get_current_status(pvgis_lookup: dict | None = None,
                       lat: float = SYSTEM_LAT,
                       lon: float = SYSTEM_LON,
                       at_time: str | None = None) -> dict | None:
    """
    Return the simulated 5-min reading for the active site.

    at_time=None      -> LIVE: reading closest to current IST wall-clock time.
    at_time="HH:MM"   -> PRESET: reading at that time-of-day (today's curve),
                         letting the dashboard show a daytime slot regardless of
                         the real clock.
    """
    if pvgis_lookup is None:
        pvgis_lookup = _get_lookup()
    utc_now  = datetime.now(timezone.utc)
    now      = utc_now + timedelta(hours=SYSTEM_UTC_OFFSET)
    readings = _readings_for_date(now.date(), pvgis_lookup, _cfg.FAULT_MODE, lat, lon)

    if at_time:
        h, m = map(int, at_time.split(":"))
        minutes = h * 60 + m
    else:
        minutes = now.hour * 60 + now.minute

    slot = min(max(minutes // 5, 0), len(readings) - 1)
    return readings[slot]


def get_intraday(target_date: date | None = None,
                 pvgis_lookup: dict | None = None,
                 lat: float = SYSTEM_LAT,
                 lon: float = SYSTEM_LON) -> list[dict]:
    """Return all 288 x 5-min readings for target_date."""
    if target_date is None:
        target_date = date.today()
    if pvgis_lookup is None:
        pvgis_lookup = _get_lookup()
    return _readings_for_date(target_date, pvgis_lookup, _cfg.FAULT_MODE, lat, lon)


def get_history(days: int = 30,
                pvgis_lookup: dict | None = None,
                lat: float = SYSTEM_LAT,
                lon: float = SYSTEM_LON) -> list[dict]:
    """Return daily summary records for the past N days (fault-free baseline)."""
    if pvgis_lookup is None:
        pvgis_lookup = _get_lookup()
    records = []
    today   = date.today()

    for d in range(days, 0, -1):
        target   = today - timedelta(days=d)
        readings = _readings_for_date(target, pvgis_lookup, "none", lat, lon)

        total_wh   = max((r["energy_wh"] or 0) for r in readings)
        peak_power = max((r["power_w"]   or 0) for r in readings)
        efficiency = round(total_wh / (SYSTEM_PEAK_KW * 1000), 3) if total_wh else 0

        records.append({
            "date":         target.strftime("%Y%m%d"),
            "energy_wh":    round(total_wh),
            "peak_power_w": round(peak_power),
            "efficiency":   efficiency,
        })

    return records


# ── Runtime fault injection ───────────────────────────────────────────────────

def inject_fault(mode: str):
    """Switch fault mode at runtime. Called by /api/fault endpoint."""
    import config
    config.FAULT_MODE = mode