"""
Fluxara Digital Twin — Mass & Energy Balance Simulation Core
Version: 1.0 · DPR v15 · September 2026

Runs as a background service alongside SCADA.
Every 15 minutes: reads SCADA state, runs balance, reports deviations.

Install: pip install numpy influxdb-client requests python-dotenv
"""

import numpy as np
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Tuple
from datetime import datetime
import json

# ═══════════════════════════════════════════════════
# CONSTANTS — From CLAUDE.md DPR v15 (LOCKED)
# ═══════════════════════════════════════════════════

@dataclass
class FluxaraConstants:
    """Locked DPR v15 constants. Do not change without full DPR revision."""
    # Operating
    operating_days_per_year: int = 330
    operating_hours_per_day: int = 24
    shifts_per_day: int = 3

    # RHA / Rice husk
    rice_husk_tpd: float = 20.0      # MT/day
    rha_purchased_tpd: float = 11.432 # MT/day
    rha_bonus_tpd: float = 3.600     # MT/day (from combustion)
    rha_total_tpd: float = 15.032    # MT/day
    rha_sio2_content: float = 0.92   # 92% SiO2
    sio2_extraction_eff: float = 0.88 # 88% extraction design

    # Products
    ps_tpd: float = 13.149           # MT/day PS Phase 1A
    pcc_tpd: float = 20.272          # MT/day PCC Phase 1A

    # Chemistry (IUPAC 2021 MW)
    MW_SiO2: float = 60.08
    MW_NaOH: float = 40.00
    MW_Na2SiO3: float = 122.06
    MW_CO2: float = 44.01
    MW_Na2CO3: float = 105.99
    MW_CaOH2: float = 74.09
    MW_CaCO3: float = 100.09
    MW_CaO: float = 56.08

    # Consumables
    naoh_fresh_tpd: float = 4.375    # MT/day fresh NaOH
    naoh_recovery_rate: float = 0.82  # 82% NaOH recovery
    cao_tpd: float = 13.364          # MT/day CaO (v15 corrected)
    cao_purity: float = 0.85         # 85% min purity
    ctab_tpd: float = 0.289          # MT/day
    stearic_tpd: float = 0.162       # MT/day

    # Energy
    gross_thermal_gjd: float = 243.0  # GJ/day gross thermal
    boiler_efficiency: float = 0.82   # 82%
    useful_heat_gjd: float = 199.3   # GJ/day after boiler
    process_heat_gjd: float = 116.3  # GJ/day with MVR
    thermal_surplus_gjd: float = 83.0 # GJ/day surplus
    mvr_saving_gjd: float = 108.6    # GJ/day saved by MVR

    # Financial (Phase 1A, 100% utilisation)
    ps_price_standard: float = 26000  # ₹/MT
    ps_price_hds: float = 47500      # ₹/MT (midpoint)
    ps_price_dental: float = 92500   # ₹/MT (midpoint)
    pcc_price_coatings: float = 11500 # ₹/MT (midpoint)
    pcc_price_sealant: float = 26000  # ₹/MT
    naoh_price_lye: float = 18240    # ₹/MT as 48% lye
    cao_price: float = 6000          # ₹/MT

CONST = FluxaraConstants()


# ═══════════════════════════════════════════════════
# PLANT STATE — what SCADA reads every 30 seconds
# ═══════════════════════════════════════════════════

@dataclass
class ScadaReading:
    """
    Real-time SCADA reading snapshot.
    All fields populated from SCADA OPC-UA tags every 15 minutes.
    None = sensor offline or stale data.
    """
    timestamp: datetime = field(default_factory=datetime.now)
    utilisation_pct: float = 100.0  # % of design capacity currently running

    # Furnace / Boiler
    furnace_temp_C: Optional[float] = None
    husk_feed_tph: Optional[float] = None    # MT/hr
    steam_pressure_bar: Optional[float] = None
    steam_temp_C: Optional[float] = None
    rha_output_tph: Optional[float] = None
    cems_co_ppm: Optional[float] = None
    cems_pm_mgm3: Optional[float] = None

    # Leach reactors
    leach_temp_C: Optional[float] = None
    rha_feed_tph: Optional[float] = None  # all reactors combined
    naoh_feed_tph: Optional[float] = None
    silicate_concentration_gL: Optional[float] = None

    # Precipitation
    precip_pH: Optional[float] = None
    precip_temp_C: Optional[float] = None
    co2_feed_kgh: Optional[float] = None   # kg/hr CO2

    # Causticisation
    caust_temp_C: Optional[float] = None
    cao_feed_tph: Optional[float] = None
    naoh_recovered_tph: Optional[float] = None

    # Evaporator / MVR
    mvr_running: bool = True
    evap_temp_eff1_C: Optional[float] = None
    evap_product_density_gml: Optional[float] = None

    # Products (daily cumulative, reset each shift)
    ps_produced_kg: float = 0.0
    pcc_produced_kg: float = 0.0


# ═══════════════════════════════════════════════════
# DIGITAL TWIN — Mass Balance Simulation
# ═══════════════════════════════════════════════════

class FluxaraDigitalTwin:
    """
    Mass and energy balance simulation.
    Call run_balance() with current SCADA readings to get:
    - Predicted flows at each stage
    - Deviations from expected
    - Alert severity for each deviation
    - Economic impact of deviations
    """

    def __init__(self, phase: str = '1A'):
        self.phase = phase
        self.const = CONST
        self.alert_history = []

    def run_balance(self, scada: ScadaReading) -> Dict:
        """
        Run full mass and energy balance given current SCADA state.
        Returns predicted values, deviations, and alerts.
        """
        u = scada.utilisation_pct / 100.0  # utilisation factor
        results = {
            'timestamp': scada.timestamp.isoformat(),
            'utilisation_pct': scada.utilisation_pct,
            'mass_balance': self._mass_balance(u),
            'energy_balance': self._energy_balance(u, scada.mvr_running),
            'deviations': self._detect_deviations(scada, u),
            'kpis': self._compute_kpis(scada, u),
            'revenue_forecast': self._revenue_forecast(u)
        }
        results['alerts'] = self._generate_alerts(results)
        return results

    def _mass_balance(self, u: float) -> Dict:
        c = self.const
        sio2_available = c.rha_total_tpd * u * c.rha_sio2_content
        sio2_extracted = sio2_available * c.sio2_extraction_eff
        na2co3_produced = sio2_extracted * (2 * c.MW_NaOH / c.MW_SiO2) * (c.MW_Na2CO3 / (2 * c.MW_NaOH))
        pcc_theoretical = na2co3_produced * (c.MW_CaCO3 / c.MW_Na2CO3)
        naoh_recovered = na2co3_produced * (2 * c.MW_NaOH / c.MW_Na2CO3) * c.naoh_recovery_rate
        cao_required = na2co3_produced * (74.09 / 105.99) / (c.cao_purity * 74.09 / 56.08)
        co2_stoich = sio2_extracted * (c.MW_CO2 / c.MW_SiO2)
        co2_actual = co2_stoich * 1.05  # 5% excess as per DPR
        return {
            'rha_to_leach_tpd': round(c.rha_total_tpd * u, 3),
            'sio2_available_tpd': round(sio2_available, 3),
            'sio2_extracted_tpd': round(sio2_extracted, 3),
            'ps_predicted_tpd': round(sio2_extracted * (c.MW_SiO2 / c.MW_SiO2), 3),
            'pcc_predicted_tpd': round(pcc_theoretical, 3),
            'na2co3_produced_tpd': round(na2co3_produced, 3),
            'naoh_recovered_tpd': round(naoh_recovered, 3),
            'naoh_fresh_makeup_tpd': round(c.naoh_fresh_tpd * u, 3),
            'cao_required_tpd': round(cao_required, 3),
            'co2_required_tpd_stoich': round(co2_stoich, 3),
            'co2_fed_tpd_5pct_excess': round(co2_actual, 3),
            'residue_tpd': round(c.rha_total_tpd * u * (1 - c.rha_sio2_content * c.sio2_extraction_eff), 3)
        }

    def _energy_balance(self, u: float, mvr_running: bool) -> Dict:
        c = self.const
        gross_thermal = c.gross_thermal_gjd * u
        useful_heat = gross_thermal * c.boiler_efficiency
        evap_no_mvr = 167.1 * u
        mvr_saving = c.mvr_saving_gjd * u if mvr_running else 0.0
        evap_with_mvr = evap_no_mvr - mvr_saving
        process_heat_total = evap_with_mvr + (42.7 + 15.1) * u
        thermal_surplus = useful_heat - process_heat_total
        return {
            'gross_thermal_gjd': round(gross_thermal, 1),
            'useful_heat_gjd': round(useful_heat, 1),
            'mvr_saving_gjd': round(mvr_saving, 1),
            'process_heat_gjd': round(process_heat_total, 1),
            'thermal_surplus_gjd': round(thermal_surplus, 1),
            'mvr_running': mvr_running,
            'alert': 'CRITICAL_MVR_DOWN' if not mvr_running and u > 0.5 else None
        }

    def _detect_deviations(self, s: ScadaReading, u: float) -> List[Dict]:
        deviations = []

        def check(name, measured, expected, tol_pct, unit, severity='WARNING'):
            if measured is None: return
            dev_pct = (abs(measured - expected) / expected) * 100
            if dev_pct > tol_pct:
                deviations.append({
                    'parameter': name, 'measured': measured,
                    'expected': expected, 'deviation_pct': round(dev_pct, 1),
                    'unit': unit, 'severity': severity
                })

        c = self.const
        check('Furnace temperature', s.furnace_temp_C, 620, 13, '°C', 'CRITICAL')
        check('Husk feed rate', s.husk_feed_tph, c.rice_husk_tpd / 24 * u, 15, 'MT/hr')
        check('Steam pressure', s.steam_pressure_bar, 7.0, 15, 'bar')
        check('Leach temperature', s.leach_temp_C, 90, 5, '°C', 'CRITICAL')
        check('Precipitation pH', s.precip_pH, 7.0, 14, 'pH', 'CRITICAL')
        check('Causticisation temperature', s.caust_temp_C, 85, 6, '°C', 'CRITICAL')
        check('CaO feed rate', s.cao_feed_tph, c.cao_tpd / 24 * u, 10, 'MT/hr')
        check('Evaporator product density', s.evap_product_density_gml, 1.20, 5, 'g/mL')
        check('CEMS CO', s.cems_co_ppm, 100, 100, 'ppm')  # 200 ppm limit
        check('CEMS PM', s.cems_pm_mgm3, 30, 67, 'mg/Nm³')  # 50 mg/Nm³ limit
        return deviations

    def _compute_kpis(self, s: ScadaReading, u: float) -> Dict:
        c = self.const
        mb = self._mass_balance(u)
        # NaOH recovery efficiency (actual vs theoretical)
        naoh_recovery = None
        if s.naoh_recovered_tph is not None and mb['na2co3_produced_tpd'] > 0:
            theoretical_naoh_per_hr = mb['na2co3_produced_tpd'] * (2 * 40 / 105.99) / 24
            naoh_recovery = s.naoh_recovered_tph / theoretical_naoh_per_hr if theoretical_naoh_per_hr > 0 else None
        # PS yield efficiency (actual vs predicted)
        ps_eff = None
        if s.ps_produced_kg > 0:
            hrs_in_shift = 8.0
            ps_expected_kg_shift = mb['ps_predicted_tpd'] * 1000 / 3  # 3 shifts
            ps_eff = s.ps_produced_kg / ps_expected_kg_shift if ps_expected_kg_shift > 0 else None
        return {
            'naoh_recovery_efficiency': round(naoh_recovery, 4) if naoh_recovery else None,
            'ps_yield_efficiency': round(ps_eff, 4) if ps_eff else None,
            'ps_predicted_tpd': round(mb['ps_predicted_tpd'], 3),
            'pcc_predicted_tpd': round(mb['pcc_predicted_tpd'], 3),
            'cao_consumption_tpd': round(mb['cao_required_tpd'], 3),
            'thermal_surplus_gjd': self._energy_balance(u, s.mvr_running)['thermal_surplus_gjd']
        }

    def _revenue_forecast(self, u: float) -> Dict:
        c = self.const
        # Assume Year 2 grade mix: 50% HDS, 50% Standard for PS; 60% coatings, 40% sealant for PCC
        ps_tpd = c.ps_tpd * u
        pcc_tpd = c.pcc_tpd * u
        ps_rev = ps_tpd * (0.5 * c.ps_price_hds + 0.5 * c.ps_price_standard)
        pcc_rev = pcc_tpd * (0.6 * c.pcc_price_coatings + 0.4 * c.pcc_price_sealant)
        residue_rev = (c.rha_total_tpd * u * (1 - c.rha_sio2_content * c.sio2_extraction_eff)) * 900
        total_daily = (ps_rev + pcc_rev + residue_rev) / 1e5  # ₹ Lakhs/day
        return {
            'ps_revenue_lakhs_day': round(ps_rev / 1e5, 2),
            'pcc_revenue_lakhs_day': round(pcc_rev / 1e5, 2),
            'total_revenue_lakhs_day': round(total_daily, 2),
            'annualised_cr': round(total_daily * c.operating_days_per_year / 100, 2)
        }

    def _generate_alerts(self, results: Dict) -> List[Dict]:
        alerts = []
        devs = results['deviations']
        energy = results['energy_balance']

        if energy.get('alert'):
            alerts.append({'severity': 'CRITICAL', 'code': energy['alert'],
                           'message': 'MVR is DOWN at utilisation >50%. Thermal deficit imminent. Check MVR compressor immediately.',
                           'economic_impact': 'Production may need to halt within 2 hours without MVR'})

        for d in devs:
            impact = self._economic_impact(d)
            alerts.append({
                'severity': d['severity'],
                'code': d['parameter'].upper().replace(' ', '_'),
                'message': f"{d['parameter']}: measured {d['measured']:.1f}{d['unit']}, expected {d['expected']:.1f}{d['unit']} ({d['deviation_pct']:.1f}% deviation)",
                'economic_impact': impact
            })

        # Sort: CRITICAL first
        alerts.sort(key=lambda a: (0 if a['severity'] == 'CRITICAL' else 1))
        return alerts

    def _economic_impact(self, deviation: Dict) -> str:
        param = deviation['parameter']
        dev_pct = deviation['deviation_pct']
        if 'Leach temperature' in param:
            # Every 5°C below 90°C drops extraction by ~4-5%
            temp_drop = max(0, 90 - (deviation['measured'] or 90))
            eff_loss = temp_drop * 0.01  # 1% efficiency per °C below 90
            daily_rev_loss = eff_loss * CONST.ps_tpd * CONST.ps_price_hds / 1e5
            return f"~₹{daily_rev_loss:.1f}L/day lost extraction (est.)"
        elif 'pH' in param:
            return "Grade risk: precipitation pH outside HDS spec — product may be Standard grade (₹19/kg lower)"
        elif 'Causticisation' in param:
            return "NaOH recovery reduced — equivalent to ₹0.47 Cr/yr per 1% CE loss"
        elif 'CEMS CO' in param:
            return "Incomplete combustion — RHA quality degrading, TSPCB compliance risk"
        elif 'Furnace temperature' in param:
            measured = deviation.get('measured', 0)
            if measured > 700:
                return "CRITICAL: RHA likely cristobalised — batch may be unusable. CHECK immediately."
            return "Low temperature — incomplete combustion, carbon in RHA"
        return f"Deviation of {dev_pct:.1f}% from design — monitor closely"


# ═══════════════════════════════════════════════════
# TELEGRAM ALERT FORMATTER
# Formats alerts for MAIS Telegram bot notification
# ═══════════════════════════════════════════════════

class TelegramAlertFormatter:
    """Format Digital Twin alerts for Telegram MAIS bot."""

    SEVERITY_EMOJI = {'CRITICAL': '🔴', 'WARNING': '🟡', 'INFO': '🟢'}

    def format_alert_message(self, results: Dict) -> str:
        alerts = results.get('alerts', [])
        kpis = results.get('kpis', {})
        rev = results.get('revenue_forecast', {})
        ts = results.get('timestamp', '')[:16]

        if not alerts:
            return (
                f"✅ *Fluxara Plant Status — {ts}*\n"
                f"All parameters within spec.\n\n"
                f"📦 PS: {kpis.get('ps_predicted_tpd', 'N/A')} MT/day\n"
                f"🪨 PCC: {kpis.get('pcc_predicted_tpd', 'N/A')} MT/day\n"
                f"💰 Revenue: ₹{rev.get('total_revenue_lakhs_day', 'N/A')}L/day"
            )

        crits = [a for a in alerts if a['severity'] == 'CRITICAL']
        warns = [a for a in alerts if a['severity'] == 'WARNING']

        header_emoji = '🔴' if crits else '🟡'
        msg = f"{header_emoji} *Fluxara Plant Alert — {ts}*\n\n"

        if crits:
            msg += "*🚨 CRITICAL ALERTS:*\n"
            for a in crits[:3]:
                msg += f"• {a['message']}\n  ↳ {a['economic_impact']}\n"

        if warns:
            msg += "\n*⚠️ WARNINGS:*\n"
            for a in warns[:3]:
                msg += f"• {a['message']}\n"

        msg += (
            f"\n*KPIs:* PS {kpis.get('ps_predicted_tpd','?')} MT/d | "
            f"PCC {kpis.get('pcc_predicted_tpd','?')} MT/d\n"
            f"💰 Revenue: ₹{rev.get('total_revenue_lakhs_day','?')}L/day"
        )
        return msg

    def format_morning_digest(self, shift_summaries: List[Dict]) -> str:
        """Format 24-hour summary for morning digest (7am daily report)."""
        return """📊 *Fluxara Daily Digest — Morning Summary*

*Production (last 24hr):*
📦 PS produced: {ps:.1f} MT (target {ps_t:.1f})
🪨 PCC produced: {pcc:.1f} MT (target {pcc_t:.1f})
⚡ Utilisation: {util:.0f}%

*Process Health:*
🔥 Furnace: avg {furn:.0f}°C (spec 550–660°C)
⚗️ Leach: avg {leach:.0f}°C (spec 88–92°C)
📉 Causticisation: NaOH recovery est. {ce:.0f}%

*Economics:*
💰 Revenue est: ₹{rev:.1f}L
📋 Critical alerts in 24hr: {crits}

Reply 'STATUS' for live readings.""".format(
            ps=0, ps_t=CONST.ps_tpd, pcc=0, pcc_t=CONST.pcc_tpd,
            util=100, furn=620, leach=90, ce=82, rev=0, crits=0
        )


# ═══════════════════════════════════════════════════
# EXAMPLE: Run the Digital Twin on mock SCADA data
# ═══════════════════════════════════════════════════
if __name__ == '__main__':
    twin = FluxaraDigitalTwin(phase='1A')
    fmt = TelegramAlertFormatter()

    # Simulate a SCADA reading with a problem (low leach temp)
    scada = ScadaReading(
        utilisation_pct=80.0,
        furnace_temp_C=635,
        husk_feed_tph=0.67,
        steam_pressure_bar=7.2,
        leach_temp_C=82,     # LOW — should be 90°C
        precip_pH=7.1,
        caust_temp_C=88,
        cao_feed_tph=0.42,
        mvr_running=True,
        evap_product_density_gml=1.19,
        cems_co_ppm=85,
        naoh_recovered_tph=0.48
    )

    results = twin.run_balance(scada)
    msg = fmt.format_alert_message(results)

    print("=== Digital Twin Output ===")
    print(json.dumps(results['mass_balance'], indent=2))
    print("\n=== Telegram Alert Preview ===")
    print(msg)