/**
 * Asset telemetry model — the heart of the NABDH digital twin.
 *
 * `sampleAsset()` is a pure function of (asset profile, timestamp, environment,
 * perturbation). Nothing here reads a clock or a random source, which is what lets
 * NABDH claim internal consistency:
 *
 *   • recorded history  = sampleAsset(profile, t_past, …)
 *   • the live stream   = sampleAsset(profile, t_now,  …)
 *   • the Time Machine  = sampleAsset(profile, t_future, …)
 *
 * The seed script writes the same values into `sensor_readings`, so the chart line
 * coming out of the database joins the live line with no discontinuity.
 */

import type { AssetTypeKey } from '../domain/enums'
import {
  type EnvironmentSample,
  solarCapacityFactor,
  windCapacityFactor,
  localHour,
} from './environment'
import { clamp, fbm, round } from './rng'

const MS_PER_MIN = 60_000

// ─────────────────────────────────────────────────────────────────────────────
// Types
// ─────────────────────────────────────────────────────────────────────────────

export interface AssetProfile {
  code: string
  assetTypeKey: AssetTypeKey
  regionCode: string
  capacityMw: number
  voltageKv: number
  ratedTempC: number
  coolingType: string
  baselineLoadPct: number
  peakHour: number
  volatility: number
  seed: number
  healthScore: number
  agingFactor: number
  profileKey: string
  installedAt: number
  lastMaintenanceAt: number | null
  customersServed: number
  isFlexibleLoad: boolean
  hasBatterySupport: boolean
  status: string
}

/** A preventive action being applied to (or simulated on) an asset. */
export interface AppliedAction {
  type: string
  assetCode: string
  magnitudeMw: number
}

/**
 * A scenario deviation from baseline conditions. Scope is narrowed by `regionCode`
 * and/or `assetCode`; leaving both null makes it national.
 */
export interface Perturbation {
  loadDeltaPct?: number
  /**
   * Optional ramp. Real demand excursions arrive over minutes, not as a step: a process
   * starts up, a tripped feeder's load transfers across, a crowd arrives. When
   * `rampMinutes` is set, the demand-side deltas grow linearly from zero at
   * `rampStartMs` to their full value `rampMinutes` later, and the thermal lag then
   * responds to that trajectory rather than to a discontinuity.
   *
   * Ambient temperature is deliberately NOT ramped — a heat wave is a background
   * condition that is already in force, not something that arrives on a schedule.
   */
  rampStartMs?: number
  rampMinutes?: number
  solarDeltaPct?: number
  windDeltaPct?: number
  tempDeltaC?: number
  outageAssetCodes?: string[]
  batteryOutage?: boolean
  lineOutage?: boolean
  evSurgeMw?: number
  industrialMw?: number
  regionCode?: string | null
  assetCode?: string | null
  actions?: AppliedAction[]
}

export const NO_PERTURBATION: Perturbation = {}

/** How much of a ramping perturbation is in force at `ts` (1 when not ramping). */
export function perturbationIntensity(p: Perturbation, ts: number): number {
  if (!p.rampMinutes || p.rampStartMs === undefined) return 1
  if (p.rampMinutes <= 0) return 1
  return clamp((ts - p.rampStartMs) / (p.rampMinutes * MS_PER_MIN), 0, 1)
}

export interface TelemetrySample {
  ts: number
  /** Utilisation as a percentage of nameplate capacity. Negative = charging. */
  loadPct: number
  loadMw: number
  voltageKv: number
  voltageDeviationPct: number
  currentA: number
  /** Top-oil / winding temperature for transformers; body temperature otherwise. */
  tempC: number
  ambientC: number
  powerFactor: number
  activePowerMw: number
  reactivePowerMvar: number
  thdPct: number
  lineLoadingPct: number
  /** Health after applying live thermal / loading stress on top of recorded health. */
  effectiveHealth: number
  /** True when a scenario has taken this asset out of service. */
  isOutage: boolean
}

// ─────────────────────────────────────────────────────────────────────────────
// Load shape
// ─────────────────────────────────────────────────────────────────────────────

interface ShapeCoefficients {
  base: number
  a1: number
  p1: number
  a2: number
  p2: number
}

const LOAD_SHAPES: Record<string, ShapeCoefficients> = {
  urban: { base: 0.8, a1: 0.2, p1: 15, a2: 0.06, p2: 20 },
  residential: { base: 0.72, a1: 0.22, p1: 16, a2: 0.1, p2: 21 },
  commercial: { base: 0.74, a1: 0.26, p1: 13, a2: 0.04, p2: 19 },
  industrial: { base: 0.9, a1: 0.08, p1: 12, a2: 0.03, p2: 2 },
  ev: { base: 0.55, a1: 0.35, p1: 19, a2: 0.12, p2: 22 },
  transmission: { base: 0.78, a1: 0.22, p1: 15, a2: 0.05, p2: 20 },
  generation: { base: 0.82, a1: 0.18, p1: 15, a2: 0.05, p2: 20 },
}

/** Normalised demand shape for the local clock hour. */
export function loadShape(profileKey: string, hour: number, peakHour: number): number {
  const c = LOAD_SHAPES[profileKey] ?? LOAD_SHAPES.urban
  const shift = peakHour - c.p1
  const value =
    c.base +
    c.a1 * Math.cos((2 * Math.PI * (hour - c.p1 - shift)) / 24) +
    c.a2 * Math.cos((4 * Math.PI * (hour - c.p2 - shift)) / 24)
  return clamp(value, 0.32, 1.3)
}

/**
 * Cooling-driven demand sensitivity to air temperature. Above ~30 °C, air-conditioning
 * makes demand rise super-linearly — the mechanism behind most summer grid stress.
 */
export function temperatureLoadFactor(tempC: number): number {
  const excess = Math.max(0, tempC - 30)
  return 1 + 0.018 * excess + 0.0006 * excess * excess
}

/** Reference conditions the seeded `baselineLoadPct` is defined at. */
const REFERENCE_TEMP_C = 35
const REFERENCE_TEMP_FACTOR = temperatureLoadFactor(REFERENCE_TEMP_C)

/** Top-oil temperature rise at rated load, by cooling class (IEC-style figures). */
const RATED_RISE_K: Record<string, number> = {
  ONAN: 55,
  ONAF: 45,
  OFAF: 40,
  OFWF: 34,
  NONE: 22,
}

// ─────────────────────────────────────────────────────────────────────────────
// Perturbation scoping
// ─────────────────────────────────────────────────────────────────────────────

export function perturbationApplies(p: Perturbation, profile: AssetProfile): boolean {
  if (p.assetCode && p.assetCode !== profile.code) return false
  if (p.regionCode && p.regionCode !== profile.regionCode) return false
  return true
}

export function isAssetOut(p: Perturbation, profile: AssetProfile): boolean {
  if (profile.status === 'offline') return true
  if (p.outageAssetCodes?.includes(profile.code)) return true
  if (p.batteryOutage && profile.assetTypeKey === 'battery_storage' && perturbationApplies(p, profile))
    return true
  if (p.lineOutage && profile.assetTypeKey === 'transmission_line' && perturbationApplies(p, profile))
    return true
  return false
}

function actionLoadReliefMw(p: Perturbation, profile: AssetProfile): number {
  if (!p.actions?.length) return 0
  let relief = 0
  for (const action of p.actions) {
    if (action.assetCode !== profile.code) continue
    if (
      action.type === 'load_transfer' ||
      action.type === 'battery_discharge' ||
      action.type === 'demand_response' ||
      action.type === 'network_reconfigure' ||
      action.type === 'generation_redispatch'
    ) {
      relief += action.magnitudeMw
    }
  }
  return relief
}

/** Multiplier applied to the thermal rise: 1 = no relief, <1 = forced cooling engaged. */
function actionCoolingRelief(p: Perturbation, profile: AssetProfile): number {
  if (!p.actions?.length) return 1
  const boosted = p.actions.some(
    (a) => a.assetCode === profile.code && a.type === 'cooling_boost',
  )
  return boosted ? 0.82 : 1
}

// ─────────────────────────────────────────────────────────────────────────────
// Raw utilisation
// ─────────────────────────────────────────────────────────────────────────────

/**
 * Utilisation before actions, as a percentage of nameplate. Split out so the thermal
 * model can look backwards in time (transformers respond to load with a lag).
 */
function rawUtilisationPct(
  profile: AssetProfile,
  ts: number,
  env: EnvironmentSample,
  p: Perturbation,
): number {
  const scoped = perturbationApplies(p, profile)
  const hour = localHour(ts)
  const minutes = ts / MS_PER_MIN
  const drift = fbm(profile.seed, minutes / 37, 3) * profile.volatility

  switch (profile.assetTypeKey) {
    case 'solar_farm': {
      const cf = solarCapacityFactor(env)
      const delta = scoped ? (p.solarDeltaPct ?? 0) : 0
      return clamp(cf * 100 * (1 + delta / 100) * (1 + drift * 0.35), 0, 100)
    }
    case 'wind_farm': {
      const cf = windCapacityFactor(env)
      const delta = scoped ? (p.windDeltaPct ?? 0) : 0
      return clamp(cf * 100 * (1 + delta / 100) * (1 + drift * 0.6), 0, 100)
    }
    case 'battery_storage': {
      // Charge through the solar midday surplus, discharge into the evening peak.
      const solar = solarCapacityFactor(env)
      const eveningPeak = Math.exp(-((hour - 19) ** 2) / 8)
      const middayCharge = Math.exp(-((hour - 12.5) ** 2) / 10) * solar
      const dispatch = eveningPeak * 0.92 - middayCharge * 0.78
      return clamp(dispatch * 100 * (1 + drift * 0.2), -100, 100)
    }
    default: {
      const shape = loadShape(profile.profileKey, hour, profile.peakHour)
      const tempRatio = temperatureLoadFactor(env.tempC) / REFERENCE_TEMP_FACTOR
      let pct = profile.baselineLoadPct * shape * tempRatio * (1 + drift)

      if (scoped) {
        const intensity = perturbationIntensity(p, ts)
        pct *= 1 + ((p.loadDeltaPct ?? 0) * intensity) / 100
        if (p.evSurgeMw && profile.assetTypeKey === 'ev_hub') {
          pct += ((p.evSurgeMw * intensity) / Math.max(profile.capacityMw, 1)) * 100
        }
        if (p.industrialMw && profile.assetTypeKey === 'industrial_load') {
          pct += ((p.industrialMw * intensity) / Math.max(profile.capacityMw, 1)) * 100
        }
      }
      return clamp(pct, 1, 165)
    }
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Public sampler
// ─────────────────────────────────────────────────────────────────────────────

export function sampleAsset(
  profile: AssetProfile,
  ts: number,
  env: EnvironmentSample,
  p: Perturbation = NO_PERTURBATION,
): TelemetrySample {
  const outage = isAssetOut(p, profile)

  if (outage) {
    return {
      ts,
      loadPct: 0,
      loadMw: 0,
      voltageKv: 0,
      voltageDeviationPct: -100,
      currentA: 0,
      tempC: env.tempC,
      ambientC: env.tempC,
      powerFactor: 0,
      activePowerMw: 0,
      reactivePowerMvar: 0,
      thdPct: 0,
      lineLoadingPct: 0,
      effectiveHealth: profile.healthScore,
      isOutage: true,
    }
  }

  const rawPct = rawUtilisationPct(profile, ts, env, p)
  const reliefMw = actionLoadReliefMw(p, profile)
  const reliefPct = (reliefMw / Math.max(profile.capacityMw, 1)) * 100
  const loadPct = clamp(rawPct - reliefPct, -100, 170)
  const loadMw = (loadPct / 100) * profile.capacityMw
  const loadPu = Math.abs(loadPct) / 100

  // Power factor sags as loading climbs past ~85 %.
  const powerFactor = clamp(
    0.985 - 0.07 * Math.max(0, loadPu - 0.85) * 4 + fbm(profile.seed + 3, ts / MS_PER_MIN / 90, 2) * 0.006,
    0.82,
    0.999,
  )

  // Voltage droop under load, plus a small stochastic component.
  const droop = 0.035 * (loadPu - 0.6)
  const voltageNoise = fbm(profile.seed + 7, ts / MS_PER_MIN / 23, 3) * 0.004
  const voltageKv = profile.voltageKv * clamp(1 - droop + voltageNoise, 0.86, 1.08)
  const voltageDeviationPct = ((voltageKv - profile.voltageKv) / profile.voltageKv) * 100

  const currentA =
    voltageKv > 0.01
      ? (Math.abs(loadMw) * 1_000_000) / (Math.sqrt(3) * voltageKv * 1000 * powerFactor)
      : 0

  // Thermal model with lag: oil temperature follows a weighted history of loading,
  // so a spike takes tens of minutes to show up as heat — which is precisely the
  // window NABDH is trying to exploit.
  const laggedPu =
    0.55 * loadPu +
    0.3 * (Math.abs(rawUtilisationPct(profile, ts - 20 * MS_PER_MIN, env, p) - reliefPct) / 100) +
    0.15 * (Math.abs(rawUtilisationPct(profile, ts - 50 * MS_PER_MIN, env, p) - reliefPct) / 100)

  const thermalTypes: AssetTypeKey[] = ['transformer', 'substation', 'power_plant']
  const ratedRise = thermalTypes.includes(profile.assetTypeKey)
    ? (RATED_RISE_K[profile.coolingType] ?? RATED_RISE_K.ONAN)
    : RATED_RISE_K.NONE

  const agingPenalty = 1 + (1 - clamp(profile.healthScore, 0, 100) / 100) * 0.5 * profile.agingFactor
  const rise = ratedRise * Math.pow(clamp(laggedPu, 0, 1.7), 1.6) * agingPenalty
  const tempC = env.tempC + rise * actionCoolingRelief(p, profile)

  const thdPct = clamp(
    1.6 + fbm(profile.seed + 13, ts / MS_PER_MIN / 45, 2) * 0.9 + 2.4 * Math.max(0, loadPu - 0.85),
    0.4,
    12,
  )

  // Live health: recorded condition, degraded by present thermal and loading stress.
  const thermalStress = clamp((tempC - profile.ratedTempC + 15) / 30, 0, 1)
  const loadStress = clamp((loadPu - 0.85) / 0.35, 0, 1)
  const effectiveHealth = clamp(
    profile.healthScore - thermalStress * 14 - loadStress * 8,
    1,
    100,
  )

  return {
    ts,
    loadPct: round(loadPct, 2),
    loadMw: round(loadMw, 3),
    voltageKv: round(voltageKv, 3),
    voltageDeviationPct: round(voltageDeviationPct, 3),
    currentA: round(currentA, 1),
    tempC: round(tempC, 2),
    ambientC: round(env.tempC, 2),
    powerFactor: round(powerFactor, 4),
    activePowerMw: round(loadMw, 3),
    reactivePowerMvar: round(Math.abs(loadMw) * Math.tan(Math.acos(powerFactor)), 3),
    thdPct: round(thdPct, 2),
    lineLoadingPct: round(Math.abs(loadPct), 2),
    effectiveHealth: round(effectiveHealth, 1),
    isOutage: false,
  }
}

/**
 * System frequency from the national supply/demand balance. Frequency is a property of
 * the whole synchronous area, not of one asset, so it is derived once per grid snapshot.
 * Saudi Arabia operates at 60 Hz.
 */
export function systemFrequency(loadMw: number, capacityMw: number, ts: number): number {
  const ratio = capacityMw > 0 ? loadMw / capacityMw : 0
  const stress = Math.max(0, ratio - 0.9)
  const wander = fbm(9173, ts / MS_PER_MIN / 11, 2) * 0.022
  return round(clamp(60 - stress * 1.2 + wander, 58.6, 60.6), 3)
}

/** Load growth over the last hour, in percentage points — a Failure DNA input. */
export function loadGrowthRate(
  profile: AssetProfile,
  ts: number,
  envNow: EnvironmentSample,
  envThen: EnvironmentSample,
  p: Perturbation = NO_PERTURBATION,
): number {
  const now = sampleAsset(profile, ts, envNow, p).loadPct
  const before = sampleAsset(profile, ts - 60 * MS_PER_MIN, envThen, p).loadPct
  return round(now - before, 3)
}

/**
 * Remaining useful life in days, from condition and current thermal stress.
 * Transformer insulation ages roughly twice as fast for every 6 °C above its rating
 * (the Montsinger rule), which is what makes the estimate move with the weather.
 */
export function remainingUsefulLifeDays(
  profile: AssetProfile,
  sample: TelemetrySample,
): number {
  const usableHealth = Math.max(0, sample.effectiveHealth - 30)
  const overTemp = Math.max(0, sample.tempC - profile.ratedTempC)
  const agingRate = Math.pow(2, overTemp / 6) * profile.agingFactor
  // Baseline: a healthy asset loses ~0.02 health points per day.
  const dailyLoss = 0.02 * agingRate + 0.004
  return Math.max(1, Math.round(usableHealth / dailyLoss))
}
