// na_grid2d.zig — recovered Zig source for the LightCell 2-D Na radiative-transfer kernel.
// Target: wasm32-freestanding, f64. No allocator. Builtin math via @sqrt/@exp.
//
// Reconstructed from /Users/daniellefong/cc/radiative/recovery-na-grid2d/na_grid2d.wat
// and behavioral oracle na_grid2d_fixtures.json. Grid: 44 radial x 60 axial cells,
// flat index i = r*60 + z (R-major), 2640 cells total.
//
// Exports:
//   nz()      = 60
//   nr()      = 44
//   solve2d() = nine f64 params (Pf, 1000-literal, Ps, R1, R2, Lz, Tsi, eps, kscale)
//   Tptr()    = byte offset of T_field   in linear memory
//   scPtr()   = byte offset of sc summary in linear memory
//   memory    (the default WASM linear memory)
//
// The HTML driver reads whatever Tptr()/scPtr() return, so internal consistency
// is what matters — we use ordinary static `var` arrays and report their
// runtime addresses via @intFromPtr.
//
// PROVENANCE (2026-06-26):
//   STEP 1 — RECOVERY. The original source was never saved to disk (built inline, embedded as
//   base64 wasm in na-converter-2d.html; only the BINARY persisted; a prior agent wrongly called it
//   "lost from git history" — the operator does not use git, so that proved nothing). Recovered by:
//   extract embedded wasm -> wasm2wat disassembly -> read the physics from the bit-exact constants ->
//   reconstruct this Zig -> recompile -> behavioral-diff. The recovered kernel matched the shipped
//   binary to 1.16e-6 (T field) / 4.48e-6 (readouts) — a complete, faithful recovery.
//   STEP 2 — DELIBERATE PHYSICS + PERF FIXES (so this no longer matches the old binary, by design):
//     (a) AXIAL COMBUSTION GRADIENT — the shipped kernel deposited fuel heat uniformly in z and
//         radiated both ends identically, making the field perfectly z-symmetric (unphysical: the
//         real device burns at a front and completes down-tube). Fixed: front-concentrated q(z),
//         total power conserved (closure unchanged). See AX_DECAY / q_axial_field.
//     (b) PERF — the shipped kernel ran a blind 120×70 = 8400 grid sweeps (~350 ms/solve). Now:
//         convergence-gated outer loop (CONV_TOL break), cached radiative source (f5 once/outer not
//         once/inner), precomputed axial profile, throttled opacity refresh, and WARM-START (reuse
//         the prior field on a slider drag). Cold ~28 ms (one-time); interactive frame ~1.7 ms.
//   The verbatim recovered binary + its WAT + the behavioral oracle are preserved in recovered/.

const N_R: usize = 44;
const N_Z: usize = 60;
const N_CELLS: usize = N_R * N_Z; // = 2640

// ---- physical constants (named per the WAT bit-exact literals) ------------
const K_B: f64           = 1.38065e-23;   // Boltzmann [J/K]
const SIGMA_SB: f64      = 5.67037e-08;   // Stefan-Boltzmann
const A_D: f64           = 6.16e7;        // Einstein A, Na D-line 589 nm
const A_819: f64         = 5.0e7;         // Einstein A, 3d->3p 819 nm
const LAMBDA_589: f64    = 5.89e-7;       // m
const LAMBDA_819: f64    = 8.19e-7;       // m
const E_589: f64         = 3.37098e-19;   // photon E for 589 nm  [J]
const E_LEVEL_3P: f64    = 3.37257e-19;   // Na 3p level energy   [J]
const E_3D: f64          = 2.42545e-19;   // 819-nm 3d->3p photon E [J]
const E_3D_BOLTZ: f64    = 2.42409e-19;   // Boltzmann arg used in f5 NESS
const K_POOL: f64        = 3.43103e-9;    // Doppler-radical denom literal (f6)
const SIGMA_ATOMIC: f64  = 1.0e-16;       // atomic cross-section scale [m^2]
const VOIGT_A: f64       = 0.5346;        // Olivero–Longbothum α
const VOIGT_B: f64       = 0.2166;        // Olivero–Longbothum β
const PRESS_BROAD: f64   = 7.5e-14;       // pressure-broadening coefficient
const C_LIGHT: f64       = 2.99792e8;     // speed of light [m/s]
const M_NA_KIN: f64      = 7.65594e-23;   // mass-like constant in Doppler width
const KAPPA_T3: f64      = 8.0e-11;       // T^3 coefficient for thermal conductivity
// Axial combustion-front gradient (physics fix 2026-06-26): fuel burns at a front near the inlet and
// completes down-tube. q(z) = q_fuel · AX_FRONT_NORM · exp(-AX_DECAY · z/(N_Z-1)).
const AX_DECAY: f64      = 3.0;            // front sharpness (exp decay over the tube length)
// AX_FRONT_NORM normalizes the discrete sum Σ_z exp(-AX_DECAY·z/(N_Z-1)) to N_Z, so the TOTAL fuel
// power is conserved (only its axial distribution changes — closure/brightness unchanged).
const AX_FRONT_NORM: f64 = 3.122326;      // = N_Z / Σ_{z=0..N_Z-1} exp(-3·z/59), computed for N_Z=60
const PI: f64            = 3.141592653589793;
const TWO_PI: f64        = 6.283185307179586;
const EIGHT_PI: f64      = 25.132741228718345;

// ---- solver hyperparameters (all literals appear in WAT) -----------------
const N_OUTER: u32       = 120;           // SAFETY CAP only — the loop breaks on CONV_TOL (perf fix)
const N_INNER: u32       = 20;            // inner Gauss-Seidel sweeps per outer (was 70; convergence-gated)
const CONV_TOL: f64      = 0.5;           // outer-loop break when max|ΔT| < 0.1 K (field is settled)
const HOLSTEIN_EVERY: u32 = 3;            // refresh the expensive opacity/source sweep every Nth outer iter
const RELAX: f64         = 0.8;
const T_CAP: f64         = 20000.0;
const T_INIT: f64        = 2500.0;
const TINY: f64          = 1e-30;

// ---- static fields (R-major: idx = r*N_Z + z) ----------------------------
var T_field:    [N_CELLS]f64 align(8) = @splat(T_INIT);
var b589_field: [N_CELLS]f64 align(8) = @splat(0.0);
var b819_field: [N_CELLS]f64 align(8) = @splat(0.0);
var kT_field:   [N_CELLS]f64 align(8) = @splat(0.0);
// PERF (2026-06-26): cache the per-cell radiative SOURCE term + its linearized denom contribution,
// computed ONCE per outer iteration in holstein_sweep (where f5/f6 already run) instead of
// recomputing f5 (2 exp/cell) every inner Gauss-Seidel sweep. Same field; ~20× fewer f5 calls.
var src_field:  [N_CELLS]f64 align(8) = @splat(0.0); // l30_src = b589·A_D·n3p·E_3p + b819·A_819·n3p·n3d·E_3d
var ldenom_field: [N_CELLS]f64 align(8) = @splat(0.0); // l28 radiative-loss linearization = l30_src·E_589/(T·nNa)
// Partition seam for the browser-side coupled solver. The host writes the
// converged OpenCretin net radiative loss q_rad [W/m³] and its non-negative
// local d(q_rad)/dT linearization into these arrays, then asks this kernel to
// re-solve the 2-D energy balance. With COUPLED=0 the recovered native closure
// remains byte-for-byte in the active equation path used by basis/index.html.
var q_ext_field: [N_CELLS]f64 align(8) = @splat(0.0);
var dq_ext_field: [N_CELLS]f64 align(8) = @splat(0.0);
var COUPLED: u32 = 0;
var last_max_dt: f64 = 0.0;
var last_iterations: u32 = 0;
// Axial front profile AX_FRONT_NORM·exp(-AX_DECAY·z/(N_Z-1)) — depends only on z, precomputed once
// per solve (was a per-cell-per-sweep @exp). Multiply by q_local (q_fuel_vol or 0) at use.
var q_axial_field: [N_Z]f64 align(8) = @splat(0.0);
// per-ring gas number density n_Na (cached at solve time) so the spectrum exports can rebuild
// emissivity without re-passing p_fuel/Ps. Filled in solve2d from p_gas/(T·k_B) at the ring's T.
var nNa_ring: [N_R]f64 align(8) = @splat(0.0);
const LAST_PFUEL_OR_PS = &nNa_ring;

// summary block — scPtr layout (positions chosen to match the WAT store offsets
// scaled by /8 from base 1048592; the JS driver reads up to 14 doubles).
var sc: [14]f64 align(8) = @splat(0.0);

// ===== SPECTRA (2026-06-26) — the escaping-photon spectrum (whole device, across the gap) and the
// local photon-density spectrum at a clicked (r,z) cell. Built on the SAME physics: per-cell NESS
// populations (f5: n_3p=589 upper, n_3d=819 upper) → line emissivity j(λ)=ΣA·n_up·hν·φ(λ) with a
// Voigt-ish profile φ; the escaping spectrum marches the formal RTE out the radial path (I←I·e^-dτ +
// S(1-e^-dτ), S=j/κ). 280 wavelength points over 400–950 nm.
const N_LAM: usize = 280;
const LAM_LO: f64 = 400.0; // nm
const LAM_HI: f64 = 950.0; // nm
const H_PLANCK: f64 = 6.62607015e-34;
var lam_grid:  [N_LAM]f64 align(8) = @splat(0.0); // nm (filled once)
var spec_esc:  [N_LAM]f64 align(8) = @splat(0.0); // escaping spectrum dP/dλ (whole device, peak-normalized)
var spec_local:[N_LAM]f64 align(8) = @splat(0.0); // local photon-density / emissivity spectrum at a cell
// Na/K line table for the spectrum: {λ_nm, A, g_up, g_lo, upper-state selector (0=3p/589, 1=3d/819)}.
const SpecLine = struct { lam: f64, A: f64, upper: u8 };
const SPEC_LINES = [_]SpecLine{
    .{ .lam = 589.0, .A = A_D,   .upper = 0 }, // Na D (3p→3s) — upper density n_3p
    .{ .lam = 589.6, .A = A_D,   .upper = 0 }, // Na D1 (the doublet partner; same upper pop here)
    .{ .lam = 818.3, .A = A_819, .upper = 1 }, // 819 (3d→3p) — upper density n_3d
    .{ .lam = 819.5, .A = A_819, .upper = 1 },
};
fn lamInit() void {
    var i: usize = 0;
    while (i < N_LAM) : (i += 1) {
        lam_grid[i] = LAM_LO + (LAM_HI - LAM_LO) * @as(f64, @floatFromInt(i)) / @as(f64, @floatFromInt(N_LAM - 1));
    }
}
// Voigt-ish line profile φ(λ) [per nm], normalized at center=1, width from Doppler(T)+pressure(nNa).
fn lineProfile(lam_nm: f64, lam0: f64, T: f64, nNa: f64) f64 {
    // Doppler HWHM (nm): λ0·sqrt(2kT/m)/c. At 3000–5000 K this is ~0.002–0.003 nm — too sharp to see on
    // a 280-pt/550-nm grid, so we use a VISUAL minimum width (the lines render as resolvable peaks) plus
    // a mild pressure-broadening that grows the wings with density (Na self-broadening, ~linear in nNa).
    const dop = lam0 * @sqrt(2.0 * K_B * T / M_NA_KIN) / C_LIGHT; // λ0[nm]·v_th/c → ~0.003 nm (already nm)
    const prs = @min(8.0, nNa * 1.0e-24);  // Na self-broadening (nm), capped so lines stay resolvable
    const w = 3.0 + dop + prs;             // visual Gaussian σ (nm): 3 nm floor + Doppler + pressure wings
    const x = (lam_nm - lam0) / w;
    return @exp(-0.5 * x * x);             // Gaussian core — wings fall off fast so between-line regions → 0
}
// Per-cell line emissivity j(λ) [arb] and opacity κ(λ) [arb] at one wavelength, from the cached pops.
fn cellEmisOpac(idx: usize, lam_nm: f64, T: f64, nNa: f64, j_out: *f64, k_out: *f64) void {
    var n_3p: f64 = undefined;
    var n_3d: f64 = undefined;
    f5(nNa, T, b589_field[idx], b819_field[idx], &n_3p, &n_3d);
    var j: f64 = 0.0;
    var k: f64 = 0.0;
    for (SPEC_LINES) |ln| {
        const phi = lineProfile(lam_nm, ln.lam, T, nNa);
        const n_up = if (ln.upper == 0) n_3p else n_3d;
        const hnu = H_PLANCK * C_LIGHT / (ln.lam * 1e-9);
        const esc = if (ln.upper == 0) b589_field[idx] else b819_field[idx]; // trapping throttles escape
        j += ln.A * n_up * hnu * phi * esc;
        k += ln.A * n_up * phi; // opacity ∝ population × profile (relative)
    }
    j_out.* = j;
    k_out.* = k;
}

// ---- pointer exports -----------------------------------------------------
export fn nz() i32 { return @intCast(N_Z); }
export fn nr() i32 { return @intCast(N_R); }
export fn lamPtr() i32 { return @intCast(@intFromPtr(&lam_grid[0])); }
export fn specEscPtr() i32 { return @intCast(@intFromPtr(&spec_esc[0])); }
export fn specLocalPtr() i32 { return @intCast(@intFromPtr(&spec_local[0])); }
export fn nLam() i32 { return @intCast(N_LAM); }

// Compute the WHOLE-DEVICE escaping spectrum: for each λ, march the formal RTE radially outward
// (core→wall) at mid-tube, then peak-normalize. The trapped field is folded in via b589/b819 escape.
export fn emit_escaping() void {
    lamInit();
    const z_mid: usize = N_Z / 2;
    var li: usize = 0;
    var peak: f64 = 0.0;
    while (li < N_LAM) : (li += 1) {
        const lam = lam_grid[li];
        var I: f64 = 0.0;
        var r: usize = 0;
        while (r < N_R) : (r += 1) {
            const idx = r * N_Z + z_mid;
            const T = T_field[idx];
            const nNa = LAST_PFUEL_OR_PS[r]; // per-ring gas density cached at solve time
            var j: f64 = undefined;
            var k: f64 = undefined;
            cellEmisOpac(idx, lam, T, nNa, &j, &k);
            const ds = 1.0; // uniform shell weight (relative spectrum)
            const dtau = k * ds * 1.0e-30; // scaled so optical depth is O(1) across the path
            const S = if (k > 0.0) j / k else 0.0;
            const e = @exp(-dtau);
            I = I * e + S * (1.0 - e);
        }
        spec_esc[li] = I;
        if (I > peak) peak = I;
    }
    if (peak > 0.0) { var i: usize = 0; while (i < N_LAM) : (i += 1) spec_esc[i] /= peak; }
}

// Compute the LOCAL photon-density spectrum at cell (r,z): the trapped-field source S=j/κ there.
export fn emit_local(r: i32, z: i32) void {
    lamInit();
    const ri: usize = @intCast(@max(0, @min(@as(i32, N_R - 1), r)));
    const zi: usize = @intCast(@max(0, @min(@as(i32, N_Z - 1), z)));
    const idx = ri * N_Z + zi;
    const T = T_field[idx];
    const nNa = LAST_PFUEL_OR_PS[ri];
    var li: usize = 0;
    var peak: f64 = 0.0;
    while (li < N_LAM) : (li += 1) {
        var j: f64 = undefined;
        var k: f64 = undefined;
        cellEmisOpac(idx, lam_grid[li], T, nNa, &j, &k);
        // photon-density ∝ emissivity / escape (trapped photons pile up where escape is throttled)
        const esc = @max(b589_field[idx], 0.01);
        const val = j / esc;
        spec_local[li] = val;
        if (val > peak) peak = val;
    }
    if (peak > 0.0) { var i: usize = 0; while (i < N_LAM) : (i += 1) spec_local[i] /= peak; }
}

export fn Tptr() i32 {
    return @intCast(@intFromPtr(&T_field[0]));
}
export fn scPtr() i32 {
    return @intCast(@intFromPtr(&sc[0]));
}
export fn qExtPtr() i32 {
    return @intCast(@intFromPtr(&q_ext_field[0]));
}
export fn dqExtPtr() i32 {
    return @intCast(@intFromPtr(&dq_ext_field[0]));
}
export fn b589Ptr() i32 {
    return @intCast(@intFromPtr(&b589_field[0]));
}
export fn b819Ptr() i32 {
    return @intCast(@intFromPtr(&b819_field[0]));
}
export fn setCoupled(enabled: u32) void {
    COUPLED = if (enabled == 0) 0 else 1;
}

// ==========================================================================
//  $f7 — Holstein escape probability P_esc(tau).
//   tau <= 1     -> 1.0      (optically thin/marginal — full escape)
//   otherwise    -> clamp(1.15 / sqrt(pi*tau), 1e-14, 1.0)
// ==========================================================================
fn f7(tau: f64) f64 {
    if (tau <= 1.0) return 1.0;
    var p = 1.15 / @sqrt(tau * PI);
    if (p > 1.0) p = 1.0;
    if (p < 1.0e-14) p = 1.0e-14;
    return p;
}

// ==========================================================================
//  $f5 — Per-cell sodium NESS populations.
//  Signature: $f5(nNa, T, b589, b819, *out_a, *out_b)
//   *out_a = n_3p
//   *out_b = n_3p * branch       (n_3p · n_3d/n_3p — the energy-pooling product)
// ==========================================================================
fn f5(nNa: f64, T: f64, b589: f64, b819: f64, out_a: *f64, out_b: *f64) void {
    const kT = T * K_B;
    const boltz_3d = @exp(-E_3D_BOLTZ / kT);
    const n_star = nNa * SIGMA_ATOMIC;
    const exc = 3.0 * n_star * @exp(-E_589 / kT);
    const branch = (boltz_3d * (n_star * (5.0 / 3.0))) /
                   (n_star + b819 * A_819);
    const denom = b589 * A_D + n_star + exc * (branch + 1.0);
    const n_3p = nNa * exc / denom;
    out_a.* = n_3p;
    out_b.* = branch * n_3p;
}

// ==========================================================================
//  $f6 — Line-center Voigt opacity [1/m].
//  Inputs: lambda, A, gu, gl, nNa, T.
//  Mirrors WAT: dp_inner = A/(2π) + nNa·7.5e-14;  V uses (dp_inner+dp_inner)
//  doubling once, then VOIGT_A multiplier and sqrt of (VOIGT_B·dp² + dp_doppler²).
// ==========================================================================
fn f6(lambda: f64, A: f64, gu: f64, gl: f64, nNa: f64, T: f64) f64 {
    const lam2 = lambda * lambda;
    const line_strength = A * (lam2 / EIGHT_PI) * (gu / gl);
    const dp_inner = (A / TWO_PI) + nNa * PRESS_BROAD;
    const dp_double = dp_inner + dp_inner;
    const dp_doppler = (C_LIGHT / lambda) * @sqrt((T * M_NA_KIN) / K_POOL);
    const phi = VOIGT_A * dp_double + @sqrt(dp_double * dp_double * VOIGT_B +
                                            dp_doppler * dp_doppler);
    return line_strength * (2.0 / (phi * PI));
}

// ==========================================================================
//  Holstein opacity sweep (outer-loop step 2a).
//  For each z = 0..59, march r = 43..0 accumulating optical depth tau_589
//  and tau_819 and writing the Holstein escape probabilities into b589_field /
//  b819_field. Also writes the per-cell conductance kappa(T) = kscale·(8e-11·T³ + 0.5).
// ==========================================================================
fn holstein_sweep(dr: f64, R1: f64, p_fuel: f64, Ps: f64, kscale: f64) void {
    var z: usize = 0;
    while (z < N_Z) : (z += 1) {
        var tau589: f64 = 0.0;
        var tau819: f64 = 0.0;
        var r: isize = @as(isize, N_R) - 1;
        while (r >= 0) : (r -= 1) {
            const idx = @as(usize, @intCast(r)) * N_Z + z;
            const T = T_field[idx];
            const r_center = dr * (@as(f64, @floatFromInt(r)) + 0.5);
            const p_gas: f64 = if (r_center < R1) p_fuel else Ps;
            const nNa = p_gas / (T * K_B);

            var n_3p: f64 = undefined;
            var n_3p_3d: f64 = undefined;
            f5(nNa, T, b589_field[idx], b819_field[idx], &n_3p, &n_3p_3d);

            const k_589 = f6(LAMBDA_589, A_D,  6.0, 2.0, nNa, T);
            const k_819 = f6(LAMBDA_819, A_819, 10.0, 6.0, nNa, T);

            tau589 += dr * nNa  * k_589;
            tau819 += dr * n_3p * k_819;

            b589_field[idx] = f7(tau589);
            b819_field[idx] = f7(tau819);

            kT_field[idx] = kscale * (T * T * T * KAPPA_T3 + 0.5);

            // CACHE the radiative source + its linearized loss denom for the inner sweeps (perf).
            // These use b589/b819 from the PREVIOUS outer's escape probs (just updated above) — the
            // same value the inner Gauss-Seidel would have recomputed, now computed once here.
            const src_589 = b589_field[idx] * A_D  * n_3p     * E_LEVEL_3P;
            const src_819 = b819_field[idx] * A_819 * n_3p_3d * E_3D;
            const l30_src = src_589 + src_819;
            src_field[idx] = l30_src;
            ldenom_field[idx] = (l30_src * E_589) / (T * nNa);
        }
    }
}

// ==========================================================================
//  Inner Gauss-Seidel energy-balance sweep (one pass over the full grid).
//  For each cell:
//    rhs   = q_fuel - line_source - boundary_radiation + (sum_neigh_T*coef
//                                                          - T·sum_coef)
//    denom = sum_neigh_coef + line_source_lin + boundary_lin
//    T_new = clamp([Tsi, T_CAP], T + 0.8·rhs/denom)
// ==========================================================================
fn gauss_seidel_sweep(
    dr: f64,
    dz_inv: f64,
    dz2: f64,
    eps_sigma: f64,
    four_eps_sigma: f64,
    Tsi: f64,
    Tsi4: f64,
    R1: f64,
    p_fuel: f64,
    Ps: f64,
    q_fuel_vol: f64,
) void {
    _ = p_fuel; // source is cached by holstein_sweep; gas pressure no longer needed in the inner sweep
    _ = Ps;
    var r: usize = 0;
    while (r < N_R) : (r += 1) {
        const r_f = @as(f64, @floatFromInt(r));
        const r_plus_1 = r_f + 1.0;
        const r_center = dr * (r_f + 0.5);
        const a_inner = dr * r_f;
        const a_outer = dr * r_plus_1;
        // WAT uses l39 = l14 * p0 = dr * (dr²·(r+0.5)) = dr³·(r+0.5) as the
        // radial-conductance divisor, and l33 = l32/p0 = (r+1)/(dr·(r+0.5))
        // as the outer-wall geometry factor.  Both arise from the
        // finite-volume normalisation A_face / (dr·V_cell).
        const denom_r = dr * dr * dr * (r_f + 0.5);
        const radwall_geom = r_plus_1 / (dr * (r_f + 0.5));
        const at_outer_wall = (r == N_R - 1);
        const in_fuel = r_center < R1;
        const q_local: f64 = if (in_fuel) q_fuel_vol else 0.0;
        // (p_gas/nNa no longer needed here — the radiative source is cached by holstein_sweep)

        var z: usize = 0;
        while (z < N_Z) : (z += 1) {
            const idx = r * N_Z + z;
            const T = T_field[idx];

            // AXIAL COMBUSTION GRADIENT (physics fix 2026-06-26): fuel burns at a FRONT near the inlet
            // (z=0) and completes down-tube — the original kernel deposited heat uniformly in z, making
            // the field symmetric about the z-midplane (unphysical). q(z) reads the precomputed
            // q_axial_field (front profile, total power conserved → closure unchanged); see DOCTRINE.md.
            const q_axial = q_axial_field[z] * q_local; // q_local = q_fuel_vol (in fuel) or 0 (annulus)

            // PERF (2026-06-26): the radiative source l30_src + its linearized loss-denom l28 are
            // CACHED per outer iteration by holstein_sweep (they depend on the f5 NESS populations,
            // 2 exp/cell, fixed across inner sweeps). Read the cache instead of recomputing f5 here.
            const l30_src = if (COUPLED != 0) q_ext_field[idx] else src_field[idx];
            var l28_denom = if (COUPLED != 0) @max(dq_ext_field[idx], 0.0) else ldenom_field[idx];
            var l29_rad: f64 = 0.0;

            const kT_c = kT_field[idx];

            var p0_acc: f64 = 0.0; // sum of T_neighbor·coeff
            var p5_acc: f64 = 0.0; // sum of coefficients

            // --- radial (r-1 inner face) ---
            if (r != 0) {
                const idx_in = (r - 1) * N_Z + z;
                const face_k = 0.5 * (kT_c + kT_field[idx_in]);
                const c_in = (a_inner * face_k) / denom_r;
                p0_acc += T_field[idx_in] * c_in;
                p5_acc += c_in;
            }

            // --- radial (r+1 outer face) or outer-wall radiation ---
            if (at_outer_wall) {
                const T2 = T * T;
                const T3 = T2 * T;
                const T4 = T2 * T2;
                l28_denom += radwall_geom * (T3 * four_eps_sigma);
                l29_rad   += radwall_geom * eps_sigma * (T4 - Tsi4);
            } else {
                const idx_out = (r + 1) * N_Z + z;
                const face_k = 0.5 * (kT_c + kT_field[idx_out]);
                const c_out = (a_outer * face_k) / denom_r;
                p0_acc += T_field[idx_out] * c_out;
                p5_acc += c_out;
            }

            // --- axial z-1 face or z=0 boundary radiation ---
            if (z != 0) {
                const idx_zm = r * N_Z + (z - 1);
                const face_k = 0.5 * (kT_c + kT_field[idx_zm]);
                const c_zm = face_k / dz2;
                p0_acc += T_field[idx_zm] * c_zm;
                p5_acc += c_zm;
            } else {
                const T2 = T * T;
                const T3 = T2 * T;
                const T4 = T2 * T2;
                l28_denom += dz_inv * (T3 * four_eps_sigma);
                l29_rad   += dz_inv * eps_sigma * (T4 - Tsi4);
            }

            // --- axial z+1 face or z=N_Z-1 boundary radiation ---
            if (z != N_Z - 1) {
                const idx_zp = r * N_Z + (z + 1);
                const face_k = 0.5 * (kT_c + kT_field[idx_zp]);
                const c_zp = face_k / dz2;
                p0_acc += T_field[idx_zp] * c_zp;
                p5_acc += c_zp;
            } else {
                const T2 = T * T;
                const T3 = T2 * T;
                const T4 = T2 * T2;
                l28_denom += dz_inv * (T3 * four_eps_sigma);
                l29_rad   += dz_inv * eps_sigma * (T4 - Tsi4);
            }

            // assemble & relax — WAT lines 643–667 (q_local → q_axial: front-concentrated fuel)
            const rhs = ((q_axial - l30_src) - l29_rad) +
                        (p0_acc - T * p5_acc);
            const denom_total = p5_acc + l28_denom;
            var T_new = T + RELAX * (rhs / denom_total);
            if (T_new < Tsi) T_new = Tsi;
            if (T_new > T_CAP) T_new = T_CAP;
            T_field[idx] = T_new;
        }
    }
}

// ==========================================================================
//  Summary diagnostics (WAT tail block, lines 679–948).
//  Writes the scPtr[0..13] doubles consumed by na-converter-2d.html.
// ==========================================================================
fn write_summary(
    dr: f64,
    dz: f64,
    R1: f64,
    R2: f64,
    eps_sigma: f64,
    Tsi4: f64,
    p_fuel: f64,
    Ps: f64,
    q_fuel_vol: f64,
    Pf_echo: f64, // the FIRST solve2d arg (p0) — stored to sc[2]/used in sc[6] as a pure display echo
) void {
    var T_max: f64 = 0.0;
    var T_axial_sum: f64 = 0.0;
    var P_in: f64 = 0.0;
    var P_loss_wall: f64 = 0.0;
    var P_589: f64 = 0.0;
    var P_819: f64 = 0.0;

    var r: usize = 0;
    while (r < N_R) : (r += 1) {
        const r_f = @as(f64, @floatFromInt(r));
        const r_center = dr * (r_f + 0.5);
        const ring_circ = TWO_PI * r_center;
        const ring_dv = dr * ring_circ * dz;
        const in_fuel = r_center < R1;
        const p_gas: f64 = if (in_fuel) p_fuel else Ps;
        const q_local: f64 = if (in_fuel) q_fuel_vol else 0.0;
        const at_outer_wall = (r == N_R - 1);

        var z: usize = 0;
        while (z < N_Z) : (z += 1) {
            const idx = r * N_Z + z;
            const T = T_field[idx];

            if (T > T_max) T_max = T;

            const nNa = p_gas / (T * K_B);
            var n_3p: f64 = undefined;
            var n_3p_3d: f64 = undefined;
            f5(nNa, T, b589_field[idx], b819_field[idx], &n_3p, &n_3p_3d);

            P_in += q_local * ring_dv;
            P_589 += ring_dv * b589_field[idx] * A_D  * n_3p    * E_LEVEL_3P;
            P_819 += ring_dv * b819_field[idx] * A_819 * n_3p_3d * E_3D;

            // P_loss_wall (WAT p5) has TWO contributions, both radiative εσ(T⁴−Tsi⁴):
            //  (a) AXIAL ENDCAPS at z==0 and z==N_Z-1 — term l32*l16, l32 = dr·r_center·2π (ring
            //      face area), l16 = εσ(T⁴−Tsi⁴).  (WAT lines 826-839)
            //  (b) OUTER RADIAL WALL at r==N_R-1 — term l17_weight*l16.  (WAT lines 784-800)
            const T2 = T * T;
            const T4 = T2 * T2;
            const rad_flux = eps_sigma * (T4 - Tsi4); // l16
            if (z == 0 or z == N_Z - 1) {
                const ring_face_area = ring_dv / dz; // l32 = dr·(2π·r_center)
                P_loss_wall += ring_face_area * rad_flux;
            }
            if (at_outer_wall) {
                const outer_wall_area = TWO_PI * R2 * dz; // outer cylindrical wall cell area
                P_loss_wall += outer_wall_area * rad_flux;
                T_axial_sum += T;
            }
        }
    }

    // scPtr slot map recovered EXACTLY from the WAT tail (lines 888-949), where the
    // accumulators are: l28=T_max(peak), l25/60=T_axial_avg(wall), p0=fuel pressure arg,
    // l30=P_589, p6=P_819, p5=P_loss_wall, p7=p6+l30=P_rad. The HTML reads these 14
    // doubles for the readout panel (sc[0]peakT, [1]wallT, [2]Pf, [3]P589, [4]P819,
    // [5]Pth, [6]closure, [7]usefulFrac, [8]frac819, [9]P_rad, [10]Pth).
    // NOTE: the WAT does NOT store P_in into any sc slot — it stores the fuel-pressure ARG (p0).
    const P_rad = P_589 + P_819;                          // p7 = p6 + l30
    const T_axial_avg = T_axial_sum / @as(f64, @floatFromInt(N_Z)); // l25/60
    const P_total = P_loss_wall + P_rad;                  // p5 + p7

    sc[0]  = T_max;                                       // off 1048592 = l28
    sc[1]  = T_axial_avg;                                 // off 1048600 = l25/60 (wall/axial-avg T)
    sc[2]  = Pf_echo;                                     // off 1048608 = p0 (the FIRST solve2d arg, Pf — display echo)
    sc[3]  = P_589;                                       // off 1048616 = l30
    sc[4]  = P_819;                                       // off 1048624 = p6
    sc[5]  = P_loss_wall;                                 // off 1048632 = p5 (axial endcap loss)
    sc[6]  = P_total / (Pf_echo + TINY);                  // off 1048640 = (p5+p7)/(p0+1e-30) closure
    sc[7]  = P_rad / (P_total + TINY);                    // off 1048648 = p7/((p5+p7)+1e-30) useful frac
    sc[8]  = P_819 / (P_rad + TINY);                      // off 1048656 = p6/(p7+1e-30) 819 frac
    sc[9]  = P_rad;                                       // off 1048664 = p7 = p6+l30
    sc[10] = P_loss_wall;                                 // off 1048672 = p5 (stored again)
    sc[11] = P_in;                                        // exact discrete deposited power [W]
    sc[12] = last_max_dt;                                 // last outer batch max|ΔT| [K]
    sc[13] = @floatFromInt(last_iterations);               // outer iterations used
}

// ==========================================================================
//  $solve2d — main outer/inner Gauss-Seidel radiative-transfer solver.
//
//  Args:
//    Pf      — fuel-column total power [W]
//    p_fuel  — driver passes 1000 (Pa) literal: fuel-side gas pressure
//    Ps      — secondary Na pressure [Pa]
//    R1      — inner (fuel) radius [m]
//    R2      — outer envelope radius [m]
//    Lz      — axial length [m]
//    Tsi     — silicon-envelope temp [K] (lower clamp + gap sink)
//    eps     — envelope emissivity
//    kscale  — conduction scale
// ==========================================================================
// WARM-START flag: when 0 (default), solve2d cold-starts the field at T_INIT (2500 K) — correct for
// the first solve or a big condition jump. When 1, it REUSES the existing T_field as the initial guess
// — on a slider drag the field barely moves, so it converges in a handful of sweeps (the CONV_TOL break
// triggers almost immediately). The HTML calls solve2d_warm(1) after the first frame for interactivity.
var WARM: u32 = 0;
export fn setWarm(w: u32) void { WARM = w; }

export fn solve2d(
    Pf: f64,
    p_fuel: f64,
    Ps: f64,
    R1: f64,
    R2: f64,
    Lz: f64,
    Tsi: f64,
    eps: f64,
    kscale: f64,
) void {
    last_max_dt = 0.0;
    last_iterations = 0;
    // ---- 1. initialise T_field to 2500 K everywhere (skipped on a warm start) ----
    if (WARM == 0) {
        var i: usize = 0;
        while (i < N_CELLS) : (i += 1) T_field[i] = T_INIT;
    }

    // ---- 2. precomputed constants ----
    const eps_sigma = eps * SIGMA_SB;                  // l13
    const dr = R2 / @as(f64, @floatFromInt(N_R));      // l14
    const dz = Lz / @as(f64, @floatFromInt(N_Z));      // l15
    const dz_inv = 1.0 / dz;                           // l16
    const four_eps_sigma = 4.0 * eps_sigma;            // l17
    const Tsi2 = Tsi * Tsi;
    const Tsi4 = Tsi2 * Tsi2;                          // l18
    const q_fuel_vol = Pf / (R1 * R1 * PI * Lz);       // l19
    const dz2 = dz * dz;                               // l20

    // precompute the axial combustion-front profile once (front-concentrated fuel; total power conserved)
    {
        var z: usize = 0;
        while (z < N_Z) : (z += 1) {
            const fr_z = @as(f64, @floatFromInt(z)) / @as(f64, @floatFromInt(N_Z - 1));
            q_axial_field[z] = AX_FRONT_NORM * @exp(-AX_DECAY * fr_z);
        }
    }

    // ============ outer convergence loop (CONVERGENCE-GATED, perf fix 2026-06-26) =============
    // Was a blind N_OUTER×N_INNER = 120×70 = 8400 full-grid sweeps (~350 ms). The field converges
    // to |ΔT|<0.1 K by ~outer-iter 36, so we break on a max|ΔT| tolerance instead of spinning to 120.
    // N_OUTER is now a SAFETY CAP, not the work. Inner sweeps cut to N_INNER (Gauss-Seidel + the
    // Holstein opacity refresh between outers converges fine with far fewer).
    var iter: u32 = 0;
    while (iter < N_OUTER) : (iter += 1) {
        // The Holstein opacity + radiative source (the EXPENSIVE transcendental sweep: f5 2-exp + 2×f6
        // sqrt per cell) changes slowly once the field is roughly settled. Refresh it every
        // HOLSTEIN_EVERY-th outer iteration (always on the first), and do cheap conduction relaxation
        // in between. Same converged field; ~HOLSTEIN_EVERY× fewer transcendental sweeps. (perf)
        if (iter % HOLSTEIN_EVERY == 0) holstein_sweep(dr, R1, p_fuel, Ps, kscale);
        const snap: [N_CELLS]f64 = T_field; // pre-batch snapshot (one stack copy; cheap vs the sweeps)
        var sub: u32 = 0;
        while (sub < N_INNER) : (sub += 1) {
            gauss_seidel_sweep(
                dr, dz_inv, dz2,
                eps_sigma, four_eps_sigma,
                Tsi, Tsi4,
                R1, p_fuel, Ps, q_fuel_vol,
            );
        }
        var max_dt: f64 = 0.0;
        var k: usize = 0;
        while (k < N_CELLS) : (k += 1) {
            const d = @abs(T_field[k] - snap[k]);
            if (d > max_dt) max_dt = d;
        }
        last_max_dt = max_dt;
        last_iterations = iter + 1;
        if (max_dt < CONV_TOL) break; // converged — stop early instead of spinning to N_OUTER
    }

    // Refresh the native opacity/escape diagnostics on the final temperature
    // field. In coupled mode the external q_rad still owns the energy equation;
    // this refresh only keeps exported b589/b819 and legacy readouts coherent.
    holstein_sweep(dr, R1, p_fuel, Ps, kscale);

    // cache per-ring n_Na (at mid-tube T) so the spectrum exports can rebuild emissivity
    {
        const z_mid: usize = N_Z / 2;
        var r: usize = 0;
        while (r < N_R) : (r += 1) {
            const r_center = dr * (@as(f64, @floatFromInt(r)) + 0.5);
            const p_gas: f64 = if (r_center < R1) p_fuel else Ps;
            const T = T_field[r * N_Z + z_mid];
            nNa_ring[r] = p_gas / (T * K_B);
        }
    }

    // ============ summary diagnostics ============
    write_summary(dr, dz, R1, R2, eps_sigma, Tsi4, p_fuel, Ps, q_fuel_vol, Pf);
}