import { round } from './rng'

/**
 * Newton–Raphson AC load flow (5.0).
 *
 * This is the real thing: it assembles a bus admittance matrix from series impedance,
 * line charging and transformer tap ratios, then solves the non-linear power-flow
 * equations in polar form until the largest power mismatch falls under a stated
 * tolerance. It reports voltage magnitude *and angle* at every bus, real and reactive
 * flow at both ends of every branch, and the losses that fall out of the difference.
 *
 * It replaces nothing. The demand allocation in `power-flow.ts` remains the fast
 * approximation, and both are selectable — a control room wants an answer in a few
 * milliseconds while it is panning a map, and a solved one when it is about to make a
 * decision. What the two must never do is be confused for each other, so every result
 * from either carries `method`, and every result from this one also carries
 * `solverStatus`, `iterations`, `convergence` and `tolerance`.
 *
 * The Newton step is damped by a backtracking line search. An undamped step is quadratic
 * near the solution and unreliable far from it: from a flat start on a heavily loaded
 * network it can overshoot out of the basin of attraction entirely, and the mismatch that
 * was falling nicely explodes in a single iteration. Halving the step until the mismatch
 * actually improves costs a few extra evaluations and turns a diverging case into a
 * converging one.
 *
 * **A diverged solution is not a solution.** When the iteration fails to converge the
 * status says so and the caller is expected to refuse to present the numbers as a
 * result. `converged` is a boolean on the object, not an inference from the absence of
 * an error.
 *
 * Conventions, stated because power-flow code is unreadable without them:
 *
 * - Per-unit throughout, on a declared MVA base. Voltages are per unit of each bus's own
 *   nominal kV, so 1.0 pu means nominal whatever the bus.
 * - Injections are generation minus load. A load bus has negative P.
 * - Angles are radians internally and degrees on the way out.
 * - A transformer's tap is on the *from* side, which is the usual convention for the
 *   off-nominal-ratio model used below.
 */

// ─────────────────────────────────────────────────────────────────────────────
// Network model
// ─────────────────────────────────────────────────────────────────────────────

export type BusType = 'slack' | 'pv' | 'pq'

export interface Bus {
  /** Asset code. The bus is the asset; there is no parallel numbering scheme. */
  code: string
  type: BusType
  /** Nominal line-to-line voltage, used as this bus's own per-unit base. */
  baseKv: number
  /** Scheduled injection in MW. Generation positive, load negative. */
  pMw: number
  /** Scheduled reactive injection in MVAr. Ignored on PV and slack buses. */
  qMvar: number
  /** Held magnitude for PV and slack buses, per unit. Ignored on PQ. */
  vSetPu?: number
  /** Reactive limits for a PV bus, in MVAr. A PV bus that hits one converts to PQ. */
  qMinMvar?: number
  qMaxMvar?: number
  /** Shunt compensation at the bus, in MVAr at nominal voltage. Positive = capacitive. */
  shuntMvar?: number
}

export interface Branch {
  from: string
  to: string
  /** Series resistance and reactance, per unit on the system base. */
  rPu: number
  xPu: number
  /** Total line charging susceptance, per unit. Split half to each end. */
  bPu?: number
  /**
   * Off-nominal tap ratio on the `from` side. 1 means nominal; a transformer stepping
   * 132 kV down to 33 kV is modelled by the two buses' different bases, not by the tap —
   * the tap is the deviation from that nominal ratio.
   */
  tap?: number
  /** Phase shift in degrees, for a phase-shifting transformer. */
  shiftDeg?: number
  /** Thermal rating in MVA, carried through so the result can report loading. */
  ratingMva?: number
  kind?: string
  /** Out of service branches are excluded from the admittance matrix entirely. */
  outOfService?: boolean
}

export interface NetworkModel {
  buses: Bus[]
  branches: Branch[]
  /** System MVA base. 100 unless a caller has a reason. */
  baseMva: number
}

export interface SolverOptions {
  /** Largest acceptable power mismatch in per unit. 1e-6 pu is 1e-4 MVA at a 100 MVA base. */
  tolerance?: number
  maxIterations?: number
  /** Flat start unless a previous solution is handed in to warm-start from. */
  initial?: { vPu: Record<string, number>; angleDeg: Record<string, number> }
  /** Enforce generator reactive limits by switching PV buses to PQ when they bind. */
  enforceQLimits?: boolean
}

// ─────────────────────────────────────────────────────────────────────────────
// Results
// ─────────────────────────────────────────────────────────────────────────────

export type SolverStatus =
  | 'converged'
  | 'max_iterations'
  | 'diverged'
  | 'singular_jacobian'
  | 'no_slack_bus'
  | 'empty_network'

export interface BusResult {
  code: string
  type: BusType
  /** Final type after any PV→PQ conversion, so a bound generator is visible. */
  finalType: BusType
  vPu: number
  vKv: number
  angleDeg: number
  /** Solved injection. On a PQ bus this equals the schedule; on slack it is the balance. */
  pMw: number
  qMvar: number
  /** Deviation from nominal, in per cent — the figure the voltage layer colours by. */
  vDeviationPct: number
  /** True when a PV bus was clamped at a reactive limit. */
  qLimited: boolean
}

export interface BranchResult {
  from: string
  to: string
  kind: string
  /** Flow measured at the sending end, positive out of `from`. */
  pFromMw: number
  qFromMvar: number
  /** Flow at the receiving end, positive into `to`. */
  pToMw: number
  qToMvar: number
  /** Always positive; direction is carried by `sourceCode`/`sinkCode`. */
  sMva: number
  currentKa: number
  sourceCode: string
  sinkCode: string
  direction: 'forward' | 'reverse' | 'idle'
  /** Real and reactive loss on the branch — the difference between the two ends. */
  lossMw: number
  lossMvar: number
  ratingMva: number
  loadingPct: number
  tap: number
}

export interface AcPowerFlowResult {
  /** Never omitted, never anything else. A caller can switch on this alone. */
  method: 'newton_raphson'
  solverStatus: SolverStatus
  converged: boolean
  iterations: number
  /** Largest remaining power mismatch, per unit. */
  convergence: number
  tolerance: number
  baseMva: number
  /** Per-iteration mismatch, so a near-miss can be told from a wild divergence. */
  mismatchHistory: number[]
  buses: BusResult[]
  branches: BranchResult[]
  slackCode: string | null
  totalGenerationMw: number
  totalLoadMw: number
  totalLossMw: number
  totalLossMvar: number
  /** Buses outside the statutory band, worst first. */
  voltageViolations: Array<{ code: string; vPu: number; kind: 'under' | 'over' }>
  overloads: Array<{ from: string; to: string; loadingPct: number }>
  /** Wall-clock is not available to a pure function; the caller stamps duration. */
  solveNote: string
}

// ─────────────────────────────────────────────────────────────────────────────
// Complex helpers, kept local and tiny
// ─────────────────────────────────────────────────────────────────────────────

interface Complex {
  re: number
  im: number
}

function cAdd(a: Complex, b: Complex): Complex {
  return { re: a.re + b.re, im: a.im + b.im }
}

function cScale(a: Complex, k: number): Complex {
  return { re: a.re * k, im: a.im * k }
}

/** 1 / (r + jx). Returns zero admittance for a zero-impedance branch rather than NaN. */
function inverseImpedance(r: number, x: number): Complex {
  const denominator = r * r + x * x
  if (denominator < 1e-12) return { re: 0, im: 0 }
  return { re: r / denominator, im: -x / denominator }
}

// ─────────────────────────────────────────────────────────────────────────────
// Admittance matrix
// ─────────────────────────────────────────────────────────────────────────────

export interface Ybus {
  g: number[][]
  b: number[][]
  index: Map<string, number>
  order: string[]
}

/**
 * Assemble the bus admittance matrix.
 *
 * The off-nominal ratio model: for a branch from `i` to `j` with series admittance `y`,
 * line charging `b/2` at each end and tap `t` on the `i` side,
 *
 *     Y_ii += y/t² + jb/2      Y_ij -= y/t
 *     Y_ji -= y/t              Y_jj += y + jb/2
 *
 * A phase shift makes the two off-diagonal terms complex conjugates of each other rather
 * than equal, which is why they are computed separately below even though they are the
 * same number whenever `shiftDeg` is zero.
 */
export function buildYbus(model: NetworkModel): Ybus {
  const order = model.buses.map((bus) => bus.code)
  const index = new Map(order.map((code, position) => [code, position]))
  const n = order.length

  const g: number[][] = Array.from({ length: n }, () => new Array<number>(n).fill(0))
  const b: number[][] = Array.from({ length: n }, () => new Array<number>(n).fill(0))

  for (const branch of model.branches) {
    if (branch.outOfService) continue
    const i = index.get(branch.from)
    const j = index.get(branch.to)
    if (i === undefined || j === undefined || i === j) continue

    const y = inverseImpedance(branch.rPu, branch.xPu)
    if (y.re === 0 && y.im === 0) continue

    const tap = branch.tap && branch.tap > 0.01 ? branch.tap : 1
    const shift = ((branch.shiftDeg ?? 0) * Math.PI) / 180
    const charging = (branch.bPu ?? 0) / 2

    // Series terms scaled by the tap.
    const yOverT2 = cScale(y, 1 / (tap * tap))
    const magnitude = 1 / tap
    const cos = Math.cos(shift)
    const sin = Math.sin(shift)

    // -y/(t e^{jθ}) for the i→j term and its conjugate partner for j→i.
    const offIj: Complex = {
      re: -magnitude * (y.re * cos + y.im * sin),
      im: -magnitude * (y.im * cos - y.re * sin),
    }
    const offJi: Complex = {
      re: -magnitude * (y.re * cos - y.im * sin),
      im: -magnitude * (y.im * cos + y.re * sin),
    }

    g[i][i] += yOverT2.re
    b[i][i] += yOverT2.im + charging
    g[j][j] += y.re
    b[j][j] += y.im + charging

    g[i][j] += offIj.re
    b[i][j] += offIj.im
    g[j][i] += offJi.re
    b[j][i] += offJi.im
  }

  // Bus shunts. A capacitor is positive MVAr and so positive susceptance.
  for (const bus of model.buses) {
    const position = index.get(bus.code)
    if (position === undefined || !bus.shuntMvar) continue
    b[position][position] += bus.shuntMvar / model.baseMva
  }

  return { g, b, index, order }
}

// ─────────────────────────────────────────────────────────────────────────────
// Linear solve
// ─────────────────────────────────────────────────────────────────────────────

/**
 * Gaussian elimination with partial pivoting.
 *
 * Dense rather than sparse: the networks this platform models are hundreds of buses, not
 * tens of thousands, and a dense factorisation of a few hundred rows costs a handful of
 * milliseconds server-side. Returns null on a singular matrix so the caller can report
 * `singular_jacobian` instead of propagating NaN into a voltage.
 */
export function solveLinear(matrix: number[][], rhs: number[]): number[] | null {
  const n = rhs.length
  if (n === 0) return []

  // Work on copies; the Jacobian is rebuilt each iteration but the caller may keep it.
  const a = matrix.map((row) => row.slice())
  const x = rhs.slice()

  for (let column = 0; column < n; column += 1) {
    let pivotRow = column
    let best = Math.abs(a[column][column])
    for (let row = column + 1; row < n; row += 1) {
      const candidate = Math.abs(a[row][column])
      if (candidate > best) {
        best = candidate
        pivotRow = row
      }
    }
    if (best < 1e-12) return null

    if (pivotRow !== column) {
      const swap = a[column]
      a[column] = a[pivotRow]
      a[pivotRow] = swap
      const value = x[column]
      x[column] = x[pivotRow]
      x[pivotRow] = value
    }

    const pivot = a[column][column]
    for (let row = column + 1; row < n; row += 1) {
      const factor = a[row][column] / pivot
      if (factor === 0) continue
      for (let k = column; k < n; k += 1) a[row][k] -= factor * a[column][k]
      x[row] -= factor * x[column]
    }
  }

  // Back substitution.
  const solution = new Array<number>(n).fill(0)
  for (let row = n - 1; row >= 0; row -= 1) {
    let sum = x[row]
    for (let k = row + 1; k < n; k += 1) sum -= a[row][k] * solution[k]
    solution[row] = sum / a[row][row]
  }

  return solution.every((value) => Number.isFinite(value)) ? solution : null
}

// ─────────────────────────────────────────────────────────────────────────────
// The solver
// ─────────────────────────────────────────────────────────────────────────────

const DEFAULT_TOLERANCE = 1e-6
const DEFAULT_MAX_ITERATIONS = 30

/** Statutory band. Outside this a bus is reported as a violation, not silently accepted. */
const V_MIN_PU = 0.94
const V_MAX_PU = 1.06

function emptyResult(status: SolverStatus, tolerance: number, baseMva: number, note: string): AcPowerFlowResult {
  return {
    method: 'newton_raphson',
    solverStatus: status,
    converged: false,
    iterations: 0,
    convergence: Number.POSITIVE_INFINITY,
    tolerance,
    baseMva,
    mismatchHistory: [],
    buses: [],
    branches: [],
    slackCode: null,
    totalGenerationMw: 0,
    totalLoadMw: 0,
    totalLossMw: 0,
    totalLossMvar: 0,
    voltageViolations: [],
    overloads: [],
    solveNote: note,
  }
}

export function solveAcPowerFlow(
  model: NetworkModel,
  options: SolverOptions = {},
): AcPowerFlowResult {
  const tolerance = options.tolerance ?? DEFAULT_TOLERANCE
  const maxIterations = options.maxIterations ?? DEFAULT_MAX_ITERATIONS
  const baseMva = model.baseMva > 0 ? model.baseMva : 100

  if (model.buses.length === 0) {
    return emptyResult('empty_network', tolerance, baseMva, 'The network has no buses.')
  }

  const slack = model.buses.find((bus) => bus.type === 'slack')
  if (!slack) {
    // Without a reference the angles are undetermined. Picking one silently would be
    // choosing an answer on the operator's behalf.
    return emptyResult(
      'no_slack_bus',
      tolerance,
      baseMva,
      'No reference bus. Angles are undetermined without one, and the solver will not choose it.',
    )
  }

  const { g, b, index, order } = buildYbus(model)
  const n = order.length
  const busAt = model.buses

  // ── State ─────────────────────────────────────────────────────────────────
  const v = new Array<number>(n).fill(1)
  const angle = new Array<number>(n).fill(0)
  // Mutable so a PV bus that hits a reactive limit can become PQ mid-solve.
  const type: BusType[] = busAt.map((bus) => bus.type)
  const qLimited = new Array<boolean>(n).fill(false)

  for (let i = 0; i < n; i += 1) {
    const bus = busAt[i]
    if (bus.type !== 'pq' && bus.vSetPu && bus.vSetPu > 0.5) v[i] = bus.vSetPu
    if (options.initial) {
      const warmV = options.initial.vPu[bus.code]
      const warmA = options.initial.angleDeg[bus.code]
      if (Number.isFinite(warmV) && warmV > 0.5) v[i] = warmV
      if (Number.isFinite(warmA)) angle[i] = (warmA * Math.PI) / 180
    }
  }

  const pSchedule = busAt.map((bus) => bus.pMw / baseMva)
  const qSchedule = busAt.map((bus) => bus.qMvar / baseMva)

  /** Injected power at a bus for any state, per unit. */
  const injectionAt = (i: number, vState: number[], aState: number[]) => {
    let p = 0
    let q = 0
    for (let j = 0; j < n; j += 1) {
      const gij = g[i][j]
      const bij = b[i][j]
      if (gij === 0 && bij === 0) continue
      const theta = aState[i] - aState[j]
      const cos = Math.cos(theta)
      const sin = Math.sin(theta)
      p += vState[j] * (gij * cos + bij * sin)
      q += vState[j] * (gij * sin - bij * cos)
    }
    return { p: vState[i] * p, q: vState[i] * q }
  }

  /** Injected power at the current state. */
  const injection = (i: number) => injectionAt(i, v, angle)

  /**
   * The largest power mismatch a trial state leaves.
   *
   * The line search needs to compare states without committing to them, so this takes the
   * state rather than reading the solver's own.
   */
  const worstMismatchAt = (
    vState: number[],
    aState: number[],
    angleBuses: number[],
    magnitudeBuses: number[],
  ) => {
    let worst = 0
    for (const i of angleBuses) {
      worst = Math.max(worst, Math.abs(pSchedule[i] - injectionAt(i, vState, aState).p))
    }
    for (const i of magnitudeBuses) {
      worst = Math.max(worst, Math.abs(qSchedule[i] - injectionAt(i, vState, aState).q))
    }
    return worst
  }

  const mismatchHistory: number[] = []
  let iterations = 0
  let convergence = Number.POSITIVE_INFINITY
  let status: SolverStatus = 'max_iterations'

  for (let iteration = 0; iteration < maxIterations; iteration += 1) {
    iterations = iteration + 1

    // ── Reactive limit enforcement, before the mismatch is measured ─────────
    if (options.enforceQLimits) {
      for (let i = 0; i < n; i += 1) {
        if (type[i] !== 'pv') continue
        const bus = busAt[i]
        const { q } = injection(i)
        const qMvar = q * baseMva
        if (bus.qMaxMvar !== undefined && qMvar > bus.qMaxMvar) {
          type[i] = 'pq'
          qSchedule[i] = bus.qMaxMvar / baseMva
          qLimited[i] = true
        } else if (bus.qMinMvar !== undefined && qMvar < bus.qMinMvar) {
          type[i] = 'pq'
          qSchedule[i] = bus.qMinMvar / baseMva
          qLimited[i] = true
        }
      }
    }

    // ── Unknowns ────────────────────────────────────────────────────────────
    const angleBuses: number[] = []
    const magnitudeBuses: number[] = []
    for (let i = 0; i < n; i += 1) {
      if (type[i] !== 'slack') angleBuses.push(i)
      if (type[i] === 'pq') magnitudeBuses.push(i)
    }

    // ── Mismatch ────────────────────────────────────────────────────────────
    const computed = new Array<{ p: number; q: number }>(n)
    for (let i = 0; i < n; i += 1) computed[i] = injection(i)

    const mismatch: number[] = []
    for (const i of angleBuses) mismatch.push(pSchedule[i] - computed[i].p)
    for (const i of magnitudeBuses) mismatch.push(qSchedule[i] - computed[i].q)

    convergence = mismatch.reduce((worst, value) => Math.max(worst, Math.abs(value)), 0)
    mismatchHistory.push(round(convergence, 12))

    if (!Number.isFinite(convergence) || convergence > 1e6) {
      status = 'diverged'
      break
    }
    if (convergence < tolerance) {
      status = 'converged'
      break
    }

    // ── Jacobian ────────────────────────────────────────────────────────────
    const rows = angleBuses.length + magnitudeBuses.length
    if (rows === 0) {
      // Every bus is the slack bus. Nothing to solve, and the mismatch is what it is.
      status = convergence < tolerance ? 'converged' : 'max_iterations'
      break
    }

    const jacobian: number[][] = Array.from({ length: rows }, () => new Array<number>(rows).fill(0))

    for (let r = 0; r < angleBuses.length; r += 1) {
      const i = angleBuses[r]
      for (let c = 0; c < angleBuses.length; c += 1) {
        const j = angleBuses[c]
        jacobian[r][c] =
          i === j
            ? -computed[i].q - b[i][i] * v[i] * v[i]
            : v[i] * v[j] * (g[i][j] * Math.sin(angle[i] - angle[j]) - b[i][j] * Math.cos(angle[i] - angle[j]))
      }
      for (let c = 0; c < magnitudeBuses.length; c += 1) {
        const j = magnitudeBuses[c]
        jacobian[r][angleBuses.length + c] =
          i === j
            ? computed[i].p / v[i] + g[i][i] * v[i]
            : v[i] * (g[i][j] * Math.cos(angle[i] - angle[j]) + b[i][j] * Math.sin(angle[i] - angle[j]))
      }
    }

    for (let r = 0; r < magnitudeBuses.length; r += 1) {
      const i = magnitudeBuses[r]
      const rowIndex = angleBuses.length + r
      for (let c = 0; c < angleBuses.length; c += 1) {
        const j = angleBuses[c]
        jacobian[rowIndex][c] =
          i === j
            ? computed[i].p - g[i][i] * v[i] * v[i]
            : -v[i] * v[j] * (g[i][j] * Math.cos(angle[i] - angle[j]) + b[i][j] * Math.sin(angle[i] - angle[j]))
      }
      for (let c = 0; c < magnitudeBuses.length; c += 1) {
        const j = magnitudeBuses[c]
        jacobian[rowIndex][angleBuses.length + c] =
          i === j
            ? computed[i].q / v[i] - b[i][i] * v[i]
            : v[i] * (g[i][j] * Math.sin(angle[i] - angle[j]) - b[i][j] * Math.cos(angle[i] - angle[j]))
      }
    }

    const step = solveLinear(jacobian, mismatch)
    if (!step) {
      status = 'singular_jacobian'
      break
    }

    // ── Damped update ───────────────────────────────────────────────────────
    //
    // Try the full Newton step first — near the solution it is the right one and the
    // search costs a single extra evaluation. Fall back by halves only when the step
    // makes things worse, which is exactly the case an undamped solver diverges on.
    const applyStep = (alpha: number) => {
      const vTrial = v.slice()
      const aTrial = angle.slice()
      for (let r = 0; r < angleBuses.length; r += 1) aTrial[angleBuses[r]] += alpha * step[r]
      for (let r = 0; r < magnitudeBuses.length; r += 1) {
        const i = magnitudeBuses[r]
        vTrial[i] += alpha * step[angleBuses.length + r]
        // A magnitude below this is a diverging iterate, not a deep brown-out.
        if (vTrial[i] < 0.05) vTrial[i] = 0.05
      }
      return { vTrial, aTrial }
    }

    let alpha = 1
    let best = applyStep(1)
    let bestMismatch = worstMismatchAt(best.vTrial, best.aTrial, angleBuses, magnitudeBuses)

    // Four halvings is enough to rescue an overshoot without turning Newton into a crawl.
    for (let attempt = 0; attempt < 4 && bestMismatch >= convergence; attempt += 1) {
      alpha /= 2
      const trial = applyStep(alpha)
      const trialMismatch = worstMismatchAt(trial.vTrial, trial.aTrial, angleBuses, magnitudeBuses)
      if (trialMismatch < bestMismatch) {
        best = trial
        bestMismatch = trialMismatch
      }
    }

    for (let i = 0; i < n; i += 1) {
      v[i] = best.vTrial[i]
      angle[i] = best.aTrial[i]
    }
  }

  const converged = status === 'converged'

  // ── Shape the answer ──────────────────────────────────────────────────────
  const finalInjection = new Array<{ p: number; q: number }>(n)
  for (let i = 0; i < n; i += 1) finalInjection[i] = injection(i)

  const buses: BusResult[] = busAt.map((bus, i) => {
    const vPu = v[i]
    return {
      code: bus.code,
      type: bus.type,
      finalType: type[i],
      vPu: round(vPu, 6),
      vKv: round(vPu * bus.baseKv, 4),
      angleDeg: round((angle[i] * 180) / Math.PI, 6),
      // On slack and PV the solved injection is the answer; on PQ it equals the schedule
      // to within the tolerance, and reporting the solved figure keeps one code path.
      pMw: round(finalInjection[i].p * baseMva, 4),
      qMvar: round(finalInjection[i].q * baseMva, 4),
      vDeviationPct: round((vPu - 1) * 100, 4),
      qLimited: qLimited[i],
    }
  })

  const branches: BranchResult[] = []
  let totalLossMw = 0
  let totalLossMvar = 0

  for (const branch of model.branches) {
    if (branch.outOfService) continue
    const i = index.get(branch.from)
    const j = index.get(branch.to)
    if (i === undefined || j === undefined || i === j) continue

    const y = inverseImpedance(branch.rPu, branch.xPu)
    if (y.re === 0 && y.im === 0) continue

    const tap = branch.tap && branch.tap > 0.01 ? branch.tap : 1
    const shift = ((branch.shiftDeg ?? 0) * Math.PI) / 180
    const charging = (branch.bPu ?? 0) / 2

    // Complex voltages, with the tap applied on the sending side.
    const vi = { re: (v[i] / tap) * Math.cos(angle[i] - shift), im: (v[i] / tap) * Math.sin(angle[i] - shift) }
    const vj = { re: v[j] * Math.cos(angle[j]), im: v[j] * Math.sin(angle[j]) }

    // Current from i into the branch: (vi - vj) * y + vi * jb/2
    const dv = { re: vi.re - vj.re, im: vi.im - vj.im }
    const iSeries = { re: dv.re * y.re - dv.im * y.im, im: dv.re * y.im + dv.im * y.re }
    const iShuntI = { re: -vi.im * charging, im: vi.re * charging }
    const iFrom = cAdd(iSeries, iShuntI)

    // S = V * conj(I)
    const sFrom = { re: vi.re * iFrom.re + vi.im * iFrom.im, im: vi.im * iFrom.re - vi.re * iFrom.im }

    const iShuntJ = { re: -vj.im * charging, im: vj.re * charging }
    const iTo = cAdd({ re: -iSeries.re, im: -iSeries.im }, iShuntJ)
    const sTo = { re: vj.re * iTo.re + vj.im * iTo.im, im: vj.im * iTo.re - vj.re * iTo.im }

    const pFromMw = sFrom.re * baseMva
    const qFromMvar = sFrom.im * baseMva
    // Reported as positive-into-`to`, which is the sign an operator expects on a one-line.
    const pToMw = -sTo.re * baseMva
    const qToMvar = -sTo.im * baseMva

    const lossMw = pFromMw - pToMw
    const lossMvar = qFromMvar - qToMvar
    totalLossMw += lossMw
    totalLossMvar += lossMvar

    const sMva = Math.hypot(pFromMw, qFromMvar)
    const ratingMva = branch.ratingMva && branch.ratingMva > 0 ? branch.ratingMva : 0
    const sendingKv = busAt[i].baseKv * v[i]
    const currentKa = sendingKv > 0.01 ? sMva / (Math.sqrt(3) * sendingKv) : 0

    const forward = pFromMw >= 0
    const idle = Math.abs(pFromMw) < 1e-6 && Math.abs(qFromMvar) < 1e-6

    branches.push({
      from: branch.from,
      to: branch.to,
      kind: branch.kind ?? 'line',
      pFromMw: round(pFromMw, 4),
      qFromMvar: round(qFromMvar, 4),
      pToMw: round(pToMw, 4),
      qToMvar: round(qToMvar, 4),
      sMva: round(sMva, 4),
      currentKa: round(currentKa, 5),
      sourceCode: forward ? branch.from : branch.to,
      sinkCode: forward ? branch.to : branch.from,
      direction: idle ? 'idle' : forward ? 'forward' : 'reverse',
      lossMw: round(lossMw, 5),
      lossMvar: round(lossMvar, 5),
      ratingMva: round(ratingMva, 2),
      loadingPct: ratingMva > 0 ? round((sMva / ratingMva) * 100, 3) : 0,
      tap,
    })
  }

  const totalGenerationMw = round(
    buses.reduce((sum, bus) => sum + Math.max(0, bus.pMw), 0),
    3,
  )
  const totalLoadMw = round(
    buses.reduce((sum, bus) => sum + Math.max(0, -bus.pMw), 0),
    3,
  )

  return {
    method: 'newton_raphson',
    solverStatus: status,
    converged,
    iterations,
    convergence: round(convergence, 12),
    tolerance,
    baseMva,
    mismatchHistory,
    buses,
    branches,
    slackCode: slack.code,
    totalGenerationMw,
    totalLoadMw,
    totalLossMw: round(totalLossMw, 4),
    totalLossMvar: round(totalLossMvar, 4),
    // Violations are only meaningful on a solved case; reporting them from a diverged
    // state would be reporting the solver's last guess as a finding about the grid.
    voltageViolations: converged
      ? buses
          .filter((bus) => bus.vPu < V_MIN_PU || bus.vPu > V_MAX_PU)
          .map((bus) => ({
            code: bus.code,
            vPu: bus.vPu,
            kind: (bus.vPu < V_MIN_PU ? 'under' : 'over') as 'under' | 'over',
          }))
          .sort((a, b) => Math.abs(b.vPu - 1) - Math.abs(a.vPu - 1))
      : [],
    overloads: converged
      ? branches
          .filter((branch) => branch.ratingMva > 0 && branch.loadingPct > 100)
          .map((branch) => ({ from: branch.from, to: branch.to, loadingPct: branch.loadingPct }))
          .sort((a, b) => b.loadingPct - a.loadingPct)
      : [],
    solveNote: converged
      ? `Converged in ${iterations} iterations to ${convergence.toExponential(2)} pu.`
      : `Did not converge (${status}). These numbers are the solver's last iterate and must not be presented as a result.`,
  }
}
