// All shader logic is inlined into this single file.
//
// HISTORY: the shader used to be split across several naga_oil modules
// (ray_gen, geodesic, stars, disk, planets, grid, skybox) that each
// `#define_import_path singularity::...` and were pulled into this file via
// `#import singularity::...`. That compiled without error, but calling ANY
// imported function at runtime produced no fragment output — the fullscreen
// quad silently drew nothing and only the camera clear color (grey) showed.
// (Local functions worked; only cross-module imports broke.) Rather than chase
// the naga_oil composition bug, every function is inlined here. The standalone
// module files still exist on disk but are no longer imported.
#import bevy_sprite::mesh2d_vertex_output::VertexOutput

struct BlackHoleUniforms {
    eye: vec4<f32>,
    forward: vec4<f32>,
    right: vec4<f32>,
    up: vec4<f32>,
    resolution: vec2<f32>,
    time: f32,
    _pad3: f32,
    rs: f32,
    disk_inner: f32,
    disk_outer: f32,
    disk_tilt: f32,
    disk_brightness: f32,
    disk_rotation_speed: f32,
    doppler_strength: f32,
    star_intensity: f32,
    skybox_intensity: f32,
    grid_density: f32,
    doppler_enabled: u32,
    grid_enabled: u32,
    planet_count: u32,
    steps: u32,
    spin: f32,
    star_aa: u32,
    bloom_threshold: f32,
    bloom_strength: f32,
    exposure: f32,
    _pad5: f32,
    disk_half_thickness: f32,
    filament_freq: f32,
    filament_sharpness: f32,
    density_freq: f32,
    density_strength: f32,
    arm_count: f32,
    arm_tightness: f32,
    arm_strength: f32,
    disk_quality: u32,
    // Disk color mode + blackbody temp (Phase 3.2), relativistic jets.
    disk_color_mode: u32,   // 0=gradient, 1=blackbody
    disk_temp: f32,         // blackbody base temperature (Kelvin)
    jets_enabled: u32,      // 0=off, 1=on
    jets_strength: f32,     // jet brightness multiplier
    aa_samples: u32,        // per-pixel supersample count (1/2/4); >1 antialiases the higher-order lensed-image rings
    _pad6: f32,             // explicit round-up to 16-byte uniform alignment
    _pad7: f32,             //   (matches BlackHoleUniforms in material.rs byte-for-byte)
};

@group(#{MATERIAL_BIND_GROUP}) @binding(0) var<uniform> uniforms: BlackHoleUniforms;

// Per-pixel fragment coordinate, set at the top of `fragment` from
// `in.position.xy` (the @builtin(position) frag coord). Used as a deterministic
// seed for the slab sub-sample jitter and the supersample jitter so the noise
// is spatially stable (no per-frame flicker). `var<private>` is the same form
// `blur.wgsl` uses for its runtime-indexed kernel tables.
var<private> pixel_seed: vec2<f32>;

// ---------- planets storage (binding 3) ----------
struct SphereData {
    center: vec4<f32>,   // xyz = center (world space), w = radius
    color: vec4<f32>,    // xyz = color, w = emissive flag
};
@group(#{MATERIAL_BIND_GROUP}) @binding(3) var<storage, read> planets: array<SphereData>;

// ---------- optional cubemap skybox (bindings 1 & 2) ----------
@group(#{MATERIAL_BIND_GROUP}) @binding(1) var skybox: texture_cube<f32>;
@group(#{MATERIAL_BIND_GROUP}) @binding(2) var skybox_sampler: sampler;

// ====================== inlined helpers ======================

// Rotate a vector around the X axis by angle a.
fn rot_x(v: vec3<f32>, a: f32) -> vec3<f32> {
    let c = cos(a);
    let s = sin(a);
    return vec3<f32>(v.x, c * v.y - s * v.z, s * v.y + c * v.z);
}

// --- ray_gen ---
// `fov` is packed into the `.w` of `up` (Rust lays out `up: Vec3` + `fov: f32`
// as one vec4 block).
fn ray_direction(uv: vec2<f32>) -> vec3<f32> {
    let tan_half_fov = tan(uniforms.up.w * 0.5);
    let dir =
        normalize(uniforms.forward.xyz)
        + uniforms.right.xyz * (uv.x * tan_half_fov)
        + uniforms.up.xyz    * (uv.y * tan_half_fov);
    return normalize(dir);
}

// --- stars ---
fn hash13(p: vec3<f32>) -> f32 {
    var q = vec3<f32>(dot(p, vec3<f32>(127.1, 311.7, 74.7)),
                      dot(p, vec3<f32>(269.5, 183.3, 246.1)),
                      dot(p, vec3<f32>(113.5, 271.9, 124.6)));
    let h = fract(sin(q) * 43758.5453);
    return h.x;
}

fn star_color(dir: vec3<f32>, intensity: f32) -> vec3<f32> {
    // The unit-sphere direction is hashed into a 3D cell grid. A pixel's ray
    // usually sits near a cell boundary, so evaluating only the home cell clips
    // any star whose center lies in a neighbor — that clip is what made stars
    // look square even with the Gaussian AA path: the hash changes at the cell
    // edge, so the falloff could never reach zero and the visible shape was
    // just the cell's bounding box. Scan the 3×3×3 neighborhood so each star's
    // radial falloff extends across the whole cell it lives in.
    let scale = 80.0;
    let p = dir * scale;
    let base = floor(p);
    let threshold = 0.985;
    // Brightest contributor wins: `best` = (brightness, r, g, b).
    var best = vec4<f32>(0.0);
    for (var dz: i32 = -1; dz <= 1; dz = dz + 1) {
        for (var dy: i32 = -1; dy <= 1; dy = dy + 1) {
            for (var dx: i32 = -1; dx <= 1; dx = dx + 1) {
                let cell = base + vec3<f32>(f32(dx), f32(dy), f32(dz));
                let h = hash13(cell);
                if (h <= threshold) { continue; }
                let b = (h - threshold) / (1.0 - threshold);
                let center = cell + vec3<f32>(0.5);
                let dist = length(p - center); // Euclidean → round, never square
                let col = mix(vec3<f32>(0.6, 0.7, 1.0), vec3<f32>(1.0, 0.9, 0.7), b);
                // Both paths share the same brightness-driven radius and the
                // same gain, so toggling AA only softens the edge — it does not
                // change peak brightness or visible count. Before, the AA path
                // used a larger radius (up to 0.65 vs 0.5), a heavier gain
                // (4.0 vs 3.0), and a slow-decaying Gaussian whose long tail
                // lifted many dim stars above the visibility floor — that's
                // what made AA look like "bigger and more" stars.
                let radius = 0.25 + b * 0.4;
                let falloff = select(
                    smoothstep(radius, 0.0, dist),                 // hard edge
                    exp(-4.6 * dist * dist / (radius * radius)),    // soft edge
                    uniforms.star_aa != 0u);
                let bright = b * falloff;
                if (bright > best.x) {
                    best = vec4<f32>(bright, col.r, col.g, col.b);
                }
            }
        }
    }
    return best.yzw * best.x * 3.0 * intensity;
}

// --- skybox ---
// textureSampleLevel (not textureSample) is REQUIRED here: this is called from
// inside the main RK45 integration loop, whose control flow is non-uniform
// (`if (accum_alpha > 0.99) { break; }`, per-pixel early exit). The WGSL spec
// forbids textureSample outside uniform control flow because it needs
// screen-space derivatives for mip selection; Chrome's Tint enforces this and
// rejects the shader, while naga (desktop) does not — which is why this only
// crashed the web build. textureSampleLevel takes an explicit LOD (0 here: the
// skybox cubemap is sampled at full resolution, no mip minification) and is
// permitted in non-uniform flow. The visual result is identical.
fn skybox_color(dir: vec3<f32>) -> vec3<f32> {
    return textureSampleLevel(skybox, skybox_sampler, dir, 0.0).rgb;
}

// --- geodesic ---
struct Deriv {
    dpos: vec3<f32>,
    ddir: vec3<f32>,
}

fn deriv(pos: vec3<f32>, dir: vec3<f32>) -> Deriv {
    let r = length(pos);
    let rs = uniforms.rs;
    // Kerr spin. χ ∈ [0,1]; a = χ·M, M = Rs/2 = 0.5 (Rs=1).
    let chi = uniforms.spin;
    let m = 0.5;
    let a = chi * m;
    // Schwarzschild radial bending (identical to Phase 1 at χ=0).
    let h = cross(pos, dir);
    let h2 = dot(h, h);
    let r5 = max(r * r * r * r * r, 1e-6);
    let radial = -1.5 * rs * h2 / r5 * pos;
    // Frame-dragging (Lense-Thirring leading term). Spin axis = +Y.
    let spin_axis = vec3<f32>(0.0, 1.0, 0.0);
    let r3 = max(r * r * r, 1e-6);
    let drag = 2.0 * m * a / r3 * cross(spin_axis, dir);
    let accel = radial + drag;
    return Deriv(dir, accel);
}

// --- disk ---
fn disk_hit(prev: vec3<f32>, cur: vec3<f32>) -> bool {
    let y0 = prev.y;
    let y1 = cur.y;
    if (y0 * y1 > 0.0) {
        return false;
    }
    let t = y0 / (y0 - y1);
    let cross = mix(prev, cur, vec3<f32>(t));
    let r = length(vec2<f32>(cross.x, cross.z));
    return r >= uniforms.disk_inner && r <= uniforms.disk_outer;
}

// --- disk noise (domain-warped FBM) ---
fn hash33(p: vec3<f32>) -> vec3<f32> {
    let q = vec3<f32>(
        dot(p, vec3<f32>(127.1, 311.7, 74.7)),
        dot(p, vec3<f32>(269.5, 183.3, 246.1)),
        dot(p, vec3<f32>(113.5, 271.9, 124.6)),
    );
    return fract(sin(q) * 43758.5453);
}

fn value_noise3(p: vec3<f32>) -> f32 {
    let i = floor(p);
    let f = fract(p);
    let u = f * f * (3.0 - 2.0 * f);
    let n000 = hash33(i + vec3<f32>(0.0, 0.0, 0.0)).x;
    let n100 = hash33(i + vec3<f32>(1.0, 0.0, 0.0)).x;
    let n010 = hash33(i + vec3<f32>(0.0, 1.0, 0.0)).x;
    let n110 = hash33(i + vec3<f32>(1.0, 1.0, 0.0)).x;
    let n001 = hash33(i + vec3<f32>(0.0, 0.0, 1.0)).x;
    let n101 = hash33(i + vec3<f32>(1.0, 0.0, 1.0)).x;
    let n011 = hash33(i + vec3<f32>(0.0, 1.0, 1.0)).x;
    let n111 = hash33(i + vec3<f32>(1.0, 1.0, 1.0)).x;
    let nx00 = mix(n000, n100, u.x);
    let nx10 = mix(n010, n110, u.x);
    let nx01 = mix(n001, n101, u.x);
    let nx11 = mix(n011, n111, u.x);
    let nxy0 = mix(nx00, nx10, u.y);
    let nxy1 = mix(nx01, nx11, u.y);
    return mix(nxy0, nxy1, u.z);
}

fn fbm3(p: vec3<f32>, octaves: u32) -> f32 {
    // MAX_OCTAVES-with-break: the conservative WebGPU form for a runtime-
    // chosen octave count. Older drivers can miscompile non-constant loop
    // bounds; this branch is the first to pass dynamic octaves here.
    const MAX_OCTAVES = 6u;
    var sum = 0.0;
    var amp = 0.5;
    var freq = 1.0;
    for (var i: u32 = 0u; i < MAX_OCTAVES; i = i + 1u) {
        if (i >= octaves) { break; }
        sum = sum + amp * value_noise3(p * freq);
        freq = freq * 2.0;
        amp = amp * 0.5;
    }
    return sum;
}

// Ridged multifractal noise: 1 - |2n-1| turns value-noise gradients into
// sharp ridges (peak where n=0.5, zero at n=0 and n=1). Raising to
// `sharpness` thins the ridges into filaments. MAX_OCTAVES-with-break is
// the conservative WebGPU form for a runtime-chosen octave count.
fn ridged_fbm(p: vec3<f32>, octaves: u32, sharpness: f32) -> f32 {
    const MAX_OCTAVES = 6u;
    var sum = 0.0;
    var amp = 0.5;
    var freq = 1.0;
    for (var i: u32 = 0u; i < MAX_OCTAVES; i = i + 1u) {
        if (i >= octaves) { break; }
        let n = value_noise3(p * freq);
        let ridge = 1.0 - abs(2.0 * n - 1.0);
        sum = sum + amp * pow(ridge, sharpness);
        freq = freq * 2.0;
        amp = amp * 0.5;
    }
    return sum;
}

fn disk_noise(pos: vec3<f32>, t: f32) -> f32 {
    let warp = fbm3(pos * 0.8 + vec3<f32>(0.0, 0.0, t * 0.1), 3u);
    let n = fbm3(pos * 2.0 + warp * 1.5 + vec3<f32>(0.0, 0.0, t * 0.3), 4u);
    return n;
}

// Result of a disk color query: emitted radiance + opacity contribution.
// Both the volumetric and flat paths return this struct so the main loop
// can treat them uniformly.
struct DiskSample {
    color: vec3<f32>,
    density: f32,
}

// Radial temperature gradient: white-hot inner → deep-orange outer.
fn temperature_color(t: f32) -> vec3<f32> {
    return mix(vec3<f32>(1.0, 0.95, 0.85), vec3<f32>(1.0, 0.45, 0.12), clamp(t, 0.0, 1.0));
}

// Analytic blackbody color (Tanner-Helland / Mitchell Charity approximation).
// Input temp in Kelvin; output linear RGB. Ported from the others/ GLSL shader.
fn blackbody(temp: f32) -> vec3f {
    // Clamp to prevent log(0) at the event horizon (infinite redshift).
    let t = max(temp, 1.0) / 100.0;
    var r: f32;
    var g: f32;
    var b: f32;
    if (t <= 66.0) {
        r = 255.0;
        g = 99.4708025861 * log(t) - 161.1195681661;
        if (t <= 19.0) {
            b = 0.0;
        } else {
            b = 138.5177312231 * log(t - 10.0) - 305.0447927307;
        }
    } else {
        r = 329.698727446 * pow(t - 60.0, -0.1332047592);
        g = 288.1221695283 * pow(t - 60.0, -0.0755148492);
        b = 255.0;
    }
    // Formula produces sRGB; convert to linear for the HDR pipeline.
    let srgb = vec3f(r, g, b) / 255.0;
    return pow(max(srgb, vec3f(0.0)), vec3f(2.2));
}

// Radial brightness falloff (∝ 1/r² from the inner edge).
fn radial_falloff(r: f32, inner: f32) -> f32 {
    return 1.0 / pow(r / inner, 2.0);
}

// Cylindrical radius in the disk plane.
fn r_of(pos: vec3<f32>) -> f32 {
    return length(vec2<f32>(pos.x, pos.z));
}

// Relativistic Doppler beaming. `dir` is the ray direction (disk-local).
fn apply_doppler(col: vec3<f32>, pos: vec3<f32>, dir: vec3<f32>) -> vec3<f32> {
    let phi = atan2(pos.z, pos.x);
    let v_orbital = sqrt(uniforms.rs / (2.0 * r_of(pos)));
    let tangent = normalize(vec3<f32>(-sin(phi), 0.0, cos(phi)));
    let vdotn = dot(tangent * v_orbital, -dir);
    let gamma = 1.0 / sqrt(max(1.0 - v_orbital * v_orbital, 1e-4));
    if (uniforms.doppler_enabled == 0u) {
        return col;
    }
    let delta = 1.0 / (gamma * (1.0 - vdotn));
    let doppler = pow(delta, 3.0) * uniforms.doppler_strength;
    return col * doppler;
}

// Kerr equatorial 4-velocity Doppler factor δ = E_obs/E_em (Page & Thorne 1974).
// Exact: solves g_μν u^μ u^ν = -1 for circular equatorial orbits, then folds in
// the conserved photon angular momentum. Returns vec2(delta, beaming=δ^3.5).
// Used by the blackbody disk color mode; the floor guards against NaN from
// pow of a negative base.
fn kerr_doppler(pos: vec3f, dir: vec3f) -> vec2f {
    let r = r_of(pos);
    let m = 0.5;
    let a = uniforms.spin * m;
    let sqrt_m = sqrt(m);
    let sign_spin = sign(uniforms.spin + 1.0e-8);
    // 1. Keplerian angular velocity Ω = dφ/dt.
    let omega = (sign_spin * sqrt_m) / (r * sqrt(r) + a * sqrt_m);
    // 2. Equatorial metric components (θ = π/2).
    let g_tt = -(1.0 - 2.0 * m / r);
    let g_tphi = -2.0 * m * a / r;
    let g_phiphi = r * r + a * a + 2.0 * m * a * a / r;
    // 3. Time component of the circular-orbit 4-velocity.
    let u_t_sq = -(g_tt + 2.0 * omega * g_tphi + omega * omega * g_phiphi);
    let u_t = 1.0 / sqrt(max(1.0e-6, u_t_sq));
    // 4. Conserved photon angular momentum (impact-parameter mapping).
    let l_photon = pos.z * dir.x - pos.x * dir.z;
    // 5. δ = E_obs/E_em = 1 / (u_t · (1 − Ω · L_photon)).
    let delta = 1.0 / max(0.01, u_t * (1.0 - omega * l_photon));
    let beaming = max(0.01, pow(delta, 3.5));
    return vec2f(delta, beaming);
}

// Unifies the color assembly for both disk paths. Mode 0 = existing gradient +
// Newtonian apply_doppler (preserves the pre-blackbody appearance). Mode 1 =
// Novikov-Thorne radial temperature gradient × Kerr δ → blackbody color × δ^3.5
// beaming, so the approaching side shifts hotter/bluer and the receding side
// cooler/redder — the physical signature absent from the gradient mode.
fn disk_emission(r: f32, pos: vec3f, dir: vec3f, brightness: f32) -> vec3f {
    let inner = uniforms.disk_inner;
    let t = (r - inner) / (uniforms.disk_outer - inner);
    let falloff = radial_falloff(r, inner);
    if (uniforms.disk_color_mode == 0u) {
        var col = temperature_color(t) * brightness * falloff * uniforms.disk_brightness;
        return apply_doppler(col, pos, dir);
    }
    // Blackbody mode: NT radial temperature profile × Kerr Doppler shift.
    let isco_r = clamp(inner / r, 0.0, 1.0);
    let nt_factor = max(0.0, 1.0 - sqrt(isco_r));
    let radial_temp = pow(isco_r, 0.75) * pow(nt_factor, 0.25);
    let d = kerr_doppler(pos, dir);
    let temperature = uniforms.disk_temp * radial_temp * d.x;
    return blackbody(temperature) * d.y * brightness * falloff * uniforms.disk_brightness;
}

// Off-tier fallback: zero-thickness disk, single midplane sample. Density is
// driven by the same radial smoothstep falloff as the volumetric path (no flat
// 0.85 floor) so the edges feather out naturally. Returns DiskSample so the
// main loop dispatches both paths uniformly.
fn disk_color_flat(pos: vec3<f32>, dir: vec3<f32>) -> DiskSample {
    let r = r_of(pos);
    let rot = uniforms.time * uniforms.disk_rotation_speed / pow(r, 1.5);
    // Domain-warped FBM for a feathered gas texture. The Keplerian shear
    // (rot ∝ 1/r^1.5) is folded into the noise flow term so inner radii flow
    // faster than outer — correct differential rotation. Amplitude is kept mild
    // so the radial temperature gradient (inside disk_emission) dominates.
    let noise = disk_noise(vec3<f32>(pos.x * 0.3, pos.z * 0.3, rot), uniforms.time);

    let col = disk_emission(r, pos, dir, 0.8 + 0.4 * noise * noise);

    // Radial smoothstep: solid through the bulk, feathered at both edges.
    let radial = smoothstep(uniforms.disk_outer, uniforms.disk_inner, r);
    return DiskSample(vec3<f32>(col), radial * 0.9);
}

// Volumetric disk color. Models the accretion disk as a physically-structured
// density field (no flat floor): radial smoothstep falloff × Gaussian vertical
// decay × a WEAK low-frequency turbulence modulation. The radial temperature
// gradient inside disk_emission drives the dominant appearance (Gargantua's
// smooth white-hot inner → deep-orange outer look); turbulence is kept to a
// subtle ±25% texture accent so the disk reads as solid gas rather than the
// jelly-like translucent slab the old 0.55-floor + ridged-filament model gave.
//
// Noise is sampled in POLAR coordinates (r_norm, phi·freq + Keplerian flow,
// height) so what little texture there is flows tangentially — correct for a
// rotating fluid. The Keplerian flow term (`+ rot`) advects the turbulence.
fn disk_color_volumetric(pos: vec3<f32>, dir: vec3<f32>) -> DiskSample {
    let r = r_of(pos);
    let phi = atan2(pos.z, pos.x);
    let rot = uniforms.time * uniforms.disk_rotation_speed / pow(r, 1.5);

    // Polar sample coordinate: (r normalized, angle × freq + Keplerian flow,
    // height within slab). The flow term advects the noise so inner radii
    // (faster rotation) drift ahead of outer radii — differential rotation.
    // NOTE: h here is normalized by disk_half_thickness, but disk_half_thickness
    // is now an H/R RATIO (semantics changed from absolute world units), so the
    // caller scales H = r * disk_half_thickness before sampling — see the
    // per-step accumulation site in the main loop.
    let r_norm = r / uniforms.disk_inner;
    let h = pos.y / max(r * uniforms.disk_half_thickness, 1e-3);
    let sp = vec3<f32>(r_norm, phi * 2.5 + rot, h);

    // Weak low-frequency turbulence. Two octaves only (was 4–5): the goal is a
    // Gargantua-style smooth disk, so turbulence is a texture accent (±25%
    // density modulation), NOT the dominant carrier. ridged filaments and
    // logarithmic-spiral arms were removed — they produced the noisy,
    // jelly-like, striated look this model replaces.
    let turbulence = fbm3(sp * uniforms.density_freq, 2u);
    let turb_mod = mix(1.0, 0.75 + 0.5 * turbulence, 0.5);  // 0.75–1.25

    // Gaussian vertical decay: densest at the midplane, ~0 at the slab edge.
    // h is normalized to the slab half-height, so exp(-(h·h)·4) → 0.018 at the
    // rim. This replaces the old uniform-in-slab density (the cause of the
    // jelly: the whole thickness glowed at a constant floor alpha).
    let h_falloff = exp(-(h * h) * 4.0);

    // Radial smoothstep falloff (outer → 0, inner → 1): smooth inner+outer
    // edges instead of a hard annulus, matching the others/ reference disk.ts.
    let radial = smoothstep(uniforms.disk_outer, uniforms.disk_inner, r);

    let total_density = radial * h_falloff * turb_mod * uniforms.density_strength;

    // Brightness: turbulence is a mild accent so the radial temperature gradient
    // (inside disk_emission) dominates the visible luminosity. This keeps the
    // disk smooth and gradient-driven rather than streaked by ridged filaments.
    let brightness = mix(0.8, 1.2, turbulence);

    let col = disk_emission(r, pos, dir, brightness);

    return DiskSample(vec3<f32>(col), total_density);
}

// Per-step disk integration that is DECOUPLED from the RK45 step size.
//
// Why this exists: the previous "sample at new_pos, weight by step_len" form
// produced radial spokes. RK45's step size is spatially adaptive (small near
// the hole where curvature is high, large in flat regions) AND varies between
// neighbouring pixels, so weighting density by `step_len` injected a per-pixel
// brightness modulation along the radial direction → spokes, densest at the
// inner disk where steps are smallest.
//
// Fix: clip THIS step's segment to the disk slab analytically and weight each
// sample by `seg_len` — the world-space length of the ray-slab intersection,
// which depends only on the ray's incidence angle and the local slab height,
// never on the RK45 step size. Then sample N points uniformly along the
// clipped segment so the integral is well-resolved (a single midpoint sample
// under-integrates the Gaussian vertical decay and leaves the disk translucent).
//
// `prev`/`new_pos` are disk-local (caller rotates by -disk_tilt).
fn integrate_disk_segment(prev: vec3f, new_pos: vec3f, dir: vec3f,
                          accum_color: ptr<function, vec3f>,
                          accum_alpha: ptr<function, f32>) {
    // Slab height H is radius-dependent (H/R ratio). Use the step's midpoint r
    // to define H for this segment — the change in r across one RK45 step is
    // small, so this local-constant approximation is accurate enough.
    let mid = (prev + new_pos) * 0.5;
    let r_mid = r_of(mid);
    let H = r_mid * uniforms.disk_half_thickness;

    // Clip the parametric segment prev + t·(new_pos − prev), t∈[0,1], to |y|≤H.
    let dy = new_pos.y - prev.y;
    var t0 = 0.0;
    var t1 = 1.0;
    if (abs(dy) < 1e-6) {
        // Step is parallel to the disk plane: in-slab iff prev is in-slab.
        if (abs(prev.y) > H) { return; }
    } else {
        let ta = ( H - prev.y) / dy;
        let tb = (-H - prev.y) / dy;
        t0 = clamp(min(ta, tb), 0.0, 1.0);
        t1 = clamp(max(ta, tb), 0.0, 1.0);
        if (t1 <= t0) { return; }  // segment does not cross the slab
    }

    // Analytic in-slab length — the key: depends on geometry, NOT on RK45 dt.
    let seg_len = (t1 - t0) * length(new_pos - prev);

    // N samples along the clipped segment, front-to-back composite. Each sample
    // carries density × (seg_len / N) so the N samples sum to the full segment
    // integral when density is roughly constant; the Gaussian vertical decay is
    // captured by sampling at differing heights within the slab. N is a
    // compile-time constant (WGSL requires static loop bounds).
    //
    // The sub-sample position is STOCHASTIC, not on the deterministic
    // (i+0.5)/N grid. A deterministic grid aliases against the 2-D pixel grid
    // and against the per-ray RK45 step lattice → a concentric moiré on the
    // disk. Folding the pixel coordinate and the sample index into a hash and
    // using it to jitter the position within stratum i turns the moiré into
    // unstructured noise (which the HDR + bloom pipeline tolerates far better
    // than coherent bands). The jitter stays inside stratum i (amplitude 1/N),
    // so the quadrature order is preserved and the integrated density is
    // unbiased. The seed is pixel-only (no time) → the noise is spatially
    // stable, no per-frame flicker.
    const N = 4u;
    for (var i: u32 = 0u; i < N; i = i + 1u) {
        let stratum = f32(i) + hash13(vec3<f32>(pixel_seed, f32(i)));
        let t = t0 + (t1 - t0) * stratum / f32(N);
        let p = mix(prev, new_pos, vec3f(t));
        let rp = r_of(p);
        if (rp < uniforms.disk_inner || rp > uniforms.disk_outer) {
            continue;  // radial clipping: outside the annulus
        }
        let s = disk_color_volumetric(p, dir);
        let ds = s.density * seg_len / f32(N);
        *accum_color += (1.0 - *accum_alpha) * s.color * ds;
        *accum_alpha += (1.0 - *accum_alpha) * ds;
    }
}

// Relativistic jets along the spin axis (Y). Bipolar cones above/below the
// disk with 0.92c outflow beaming. Ported from others/' sample_relativistic_jets:
// Gaussian radial falloff, exponential length decay, outward-flowing noise,
// δ^3.5 beaming (no floor needed — β=0.92 keeps the denominator positive).
// Front-to-back composited into the same accumulators as the disk.
//
// `step_len` is the world-space length of THIS RK45 step (|new_pos − prev|),
// passed in by the caller. The per-step emission is weighted by it so the
// integrated jet brightness is step-size-independent — the same technique the
// disk path uses (integrate_disk_segment weights on the analytic in-slab
// length). Two prior attempts both produced banding: `* dt` varied with the
// adaptive step size (radial bands), and a fixed `0.5` made the per-UNIT-length
// contribution ∝ 1/step_len, so it brightened wherever steps packed together
// (near the hole) → the concentric base-of-jet bands. Weighting on the actual
// geometric step length makes the contribution per unit length constant.
fn sample_jets(pos: vec3f, dir: vec3f, r_plus: f32, step_len: f32,
               accum_color: ptr<function, vec3f>, accum_alpha: ptr<function, f32>) {
    // Relativistic jets are spin-powered (Blandford-Znajek): the mechanism taps
    // the ergosphere, which only exists for a rotating hole. At χ ≈ 0 there is
    // nothing to drive an outflow, so suppress the jets regardless of the
    // jets_enabled toggle — the toggle expresses user intent, spin expresses
    // physics. Without this gate the default scene (spin = 0, jets_enabled =
    // true) shows blue/white columns over the poles that have no physical cause.
    if (uniforms.spin < 0.05) { return; }
    let jet_v = abs(pos.y);
    let jet_max_h = 80.0;
    if (jet_v <= r_plus * 1.8 || jet_v >= jet_max_h) {
        return;
    }
    let jet_r = length(vec2f(pos.x, pos.z));
    let jet_width = 1.0 + jet_v * 0.15;
    if (jet_r >= jet_width * 2.0) {
        return;
    }
    let radial_falloff = exp(-(jet_r * jet_r) / (jet_width * 0.5));
    let length_falloff = exp(-jet_v * 0.05);

    // Outward-flowing turbulence: the y term reverses sign of the time flow so
    // the lower jet streams downward and the upper jet upward.
    let flow = pos.y * 2.0 - uniforms.time * 8.0;
    let uv_jet = vec3f(pos.x, flow, pos.z);
    let noise_val = value_noise3(uv_jet * 0.5) * 0.6 + value_noise3(uv_jet * 1.5) * 0.4;
    let jet_density = radial_falloff * length_falloff * max(0.0, noise_val - 0.2);

    if (jet_density <= 0.001) {
        return;
    }
    // 0.92c outflow beaming.
    let jet_vel = 0.92 * sign(pos.y);
    let jet_vel_vec = vec3f(0.0, jet_vel, 0.0);
    let cos_theta = dot(normalize(jet_vel_vec), -dir);
    let beta = abs(jet_vel);
    let gamma = 1.0 / sqrt(1.0 - beta * beta);
    let delta = 1.0 / (gamma * (1.0 - beta * cos_theta));
    // Cap the beaming: with β=0.92 the approaching jet's δ^3.5 reaches ~258×,
    // which blows the blue jet past the bloom threshold and saturates it to
    // white. Clamp to 8× — still visibly brighter on the approaching side, but
    // the base color (0.4, 0.7, 1.0) survives the HDR + bloom pipeline.
    let beaming = min(pow(delta, 3.5), 8.0);

    let base_color = vec3f(0.4, 0.7, 1.0);
    // Weight the per-step emission by the geometric step length so the
    // integrated brightness is proportional to the path length through the jet,
    // not to how many RK45 steps happened to fall there. The 0.1 coefficient is
    // a per-unit-length strength, calibrated against the default scene to match
    // the old fixed-0.5 brightness at the typical accepted-step length (~0.5).
    let emission = base_color * jet_density * 0.05 * beaming * uniforms.jets_strength * step_len * 0.1;
    *accum_color += emission * (1.0 - *accum_alpha);
    *accum_alpha += jet_density * 0.05 * uniforms.jets_strength * step_len * 0.1;
}

// --- planets ---
// `prev`/`cur` are in DISK-LOCAL space; planet centers are world space, so we
// rotate each center into disk-local space here.
fn planet_hit(prev: vec3<f32>, cur: vec3<f32>, dir: vec3<f32>) -> vec4<f32> {
    var nearest_t = 1e9;
    var nearest_col = vec3<f32>(0.0);
    var found = false;
    for (var i: u32 = 0u; i < uniforms.planet_count; i = i + 1u) {
        let s = planets[i];
        let center = rot_x(s.center.xyz, -uniforms.disk_tilt);
        let radius = s.center.w;
        let seg = cur - prev;
        let oc = prev - center;
        let a = dot(seg, seg);
        let b = 2.0 * dot(oc, seg);
        let c = dot(oc, oc) - radius * radius;
        let disc = b * b - 4.0 * a * c;
        if (disc < 0.0) { continue; }
        let sq = sqrt(disc);
        var t = (-b - sq) / (2.0 * a);
        if (t < 0.0) { t = (-b + sq) / (2.0 * a); }
        if (t >= 0.0 && t <= 1.0 && t < nearest_t) {
            nearest_t = t;
            let hit_pos = prev + seg * t;
            let n = normalize(hit_pos - center);
            let light_dir = normalize(vec3<f32>(0.5, 0.8, 0.3));
            let ndl = max(dot(n, light_dir), 0.0);
            var col = s.color.xyz * (0.2 + 0.8 * ndl);
            if (s.color.w > 0.5) { col = s.color.xyz; }
            nearest_col = col;
            found = true;
        }
    }
    if (found) {
        return vec4<f32>(nearest_col, 0.95);
    }
    return vec4<f32>(0.0, 0.0, 0.0, 0.0);
}

// --- grid (Flamm's paraboloid) ---
fn flamm_depth(r: f32) -> f32 {
    if (r <= uniforms.rs) { return 0.0; }
    return -2.0 * sqrt(uniforms.rs * (r - uniforms.rs));
}

fn grid_hit(prev: vec3<f32>, cur: vec3<f32>) -> vec3<f32> {
    let r0 = length(vec2<f32>(prev.x, prev.z));
    let r1 = length(vec2<f32>(cur.x, cur.z));
    let z0_surf = flamm_depth(r0);
    let z1_surf = flamm_depth(r1);
    if ((prev.y - z0_surf) * (cur.y - z1_surf) > 0.0) {
        return vec3<f32>(0.0);
    }
    var hit = vec3<f32>(0.0);
    var found = false;
    for (var s: i32 = 0; s < 8; s = s + 1) {
        let f = f32(s + 1) / 8.0;
        let p = mix(prev, cur, vec3<f32>(f));
        let r = length(vec2<f32>(p.x, p.z));
        let surf = flamm_depth(r);
        if (abs(p.y - surf) < 0.3) {
            hit = p;
            found = true;
            break;
        }
    }
    if (!found) { return vec3<f32>(0.0); }

    let r = length(vec2<f32>(hit.x, hit.z));
    let phi = atan2(hit.z, hit.x);
    let ring = smoothstep(0.06, 0.0, abs(fract(r * uniforms.grid_density * 0.5) - 0.5));
    let spoke = smoothstep(0.04, 0.0, abs(fract(phi * 6.0 / 6.283185) - 0.5));
    let grid = max(ring, spoke);
    let fade = smoothstep(-15.0, -1.0, hit.y);
    let col = vec3<f32>(0.15, 0.3, 0.6) * grid * fade;
    return col * 0.5;
}

// One Dormand-Prince RK45 step. Returns the 5th-order solution and the
// error estimate (y5 - y4) as a vec3 (position error; direction error is
// folded in via normalize so we only need position error for step control).
struct RkStep {
    pos: vec3<f32>,
    dir: vec3<f32>,
    err: f32,
};

fn rk45_step(pos: vec3<f32>, dir: vec3<f32>, dt: f32) -> RkStep {
    // Butcher tableau (Dormand-Prince), 6 stages. Each deriv() returns Deriv{dpos, ddir}.
    let k1 = deriv(pos, dir);
    let p2 = pos + k1.dpos * dt * 0.2;
    let d2 = normalize(dir + k1.ddir * dt * 0.2);
    let k2 = deriv(p2, d2);
    let p3 = pos + (k1.dpos * 0.075 + k2.dpos * 0.225) * dt;
    let d3 = normalize(dir + (k1.ddir * 0.075 + k2.ddir * 0.225) * dt);
    let k3 = deriv(p3, d3);
    let p4 = pos + (k1.dpos * 0.3 + k2.dpos * -0.9 + k3.dpos * 1.2) * dt;
    let d4 = normalize(dir + (k1.ddir * 0.3 + k2.ddir * -0.9 + k3.ddir * 1.2) * dt);
    let k4 = deriv(p4, d4);
    let p5 = pos + (k1.dpos * -11.0/54.0 + k2.dpos * 2.5 + k3.dpos * -70.0/27.0 + k4.dpos * 35.0/27.0) * dt;
    let d5 = normalize(dir + (k1.ddir * -11.0/54.0 + k2.ddir * 2.5 + k3.ddir * -70.0/27.0 + k4.ddir * 35.0/27.0) * dt);
    let k5 = deriv(p5, d5);
    let p6 = pos + (k1.dpos * 1631.0/55296.0 + k2.dpos * 175.0/512.0 + k3.dpos * 575.0/13824.0 + k4.dpos * 44275.0/110592.0 + k5.dpos * 253.0/4096.0) * dt;
    let d6 = normalize(dir + (k1.ddir * 1631.0/55296.0 + k2.ddir * 175.0/512.0 + k3.ddir * 575.0/13824.0 + k4.ddir * 44275.0/110592.0 + k5.ddir * 253.0/4096.0) * dt);
    let k6 = deriv(p6, d6);
    // 5th-order solution (used to advance).
    let new_pos = pos + (k1.dpos * 37.0/378.0 + k3.dpos * 250.0/621.0 + k4.dpos * 125.0/594.0 + k5.dpos * 512.0/1771.0 + k6.dpos * 0.0) * dt;
    let new_dir = normalize(dir + (k1.ddir * 37.0/378.0 + k3.ddir * 250.0/621.0 + k4.ddir * 125.0/594.0 + k5.ddir * 512.0/1771.0 + k6.ddir * 0.0) * dt);
    // 4th-order solution (for error estimate).
    let pos4 = pos + (k1.dpos * 2825.0/27648.0 + k3.dpos * 18575.0/48384.0 + k4.dpos * 13525.0/55296.0 + k5.dpos * 277.0/14336.0 + k6.dpos * 0.25) * dt;
    let err = length(new_pos - pos4);
    return RkStep(new_pos, new_dir, err);
}

// ====================== main ======================

// One full ray march (camera → capture/escape/budget) for a single ray
// direction, returning the accumulated HDR color. Factored out of `fragment`
// so the supersampling loop can fire several jittered rays per pixel and
// average them. `dir` is the camera-space ray direction; the disk-local
// rotation and all per-ray state live in here so each sample is independent.
fn march(dir: vec3<f32>) -> vec3<f32> {
    // Total path length to integrate: enough to go from the camera, past the
    // hole, and far enough beyond to count as escaped.
    let eye_dist = length(uniforms.eye.xyz);
    let escape_r = max(eye_dist * 2.0, 100.0);
    let total_path = eye_dist + escape_r;
    // Adaptive RK45 constants.
    let steps_max = uniforms.steps;
    let dt_init = total_path / f32(steps_max);
    let dt_min = dt_init * 0.25;
    let dt_max = dt_init * 4.0;
    let tol = 1e-3;
    let r_plus = 0.5 + sqrt(max(0.25 - (uniforms.spin * 0.5) * (uniforms.spin * 0.5), 0.0));

    // Work in disk-local space: rotate eye + dir by -disk_tilt around X so the
    // disk lies on y=0. (disk_hit/disk_color assume disk-local coords.)
    var pos = rot_x(uniforms.eye.xyz, -uniforms.disk_tilt);
    var d   = normalize(rot_x(dir, -uniforms.disk_tilt));
    var dt  = dt_init;
    var prev = pos;
    var budget = steps_max;

    var accum_color = vec3<f32>(0.0);
    var accum_alpha = 0.0;

    loop {
        if (budget == 0u) { break; }

        let step = rk45_step(pos, d, dt);
        let err = step.err;

        if (err > tol * 10.0 && dt > dt_min) {
            // Reject: shrink dt, retry (does not consume budget).
            // The `dt > dt_min` guard makes dt_min a forced-accept floor: once
            // dt is already at its minimum, accept the step regardless of error
            // rather than spinning forever on `continue`. This prevents a GPU
            // hang when a ray hits a region whose error is intrinsically above
            // tolerance even at the smallest step (extreme spin, near-horizon).
            dt = dt_min;
            continue;
        }
        // Accept: consume one budget unit, refine dt.
        budget = budget - 1u;
        if (err > tol * 10.0) {
            // Forced accept at dt_min (err still high): don't grow dt back up.
        } else {
            dt = clamp(dt * pow(tol / max(err, 1e-12), 0.2), dt_min, dt_max);
        }

        let new_pos = step.pos;
        let new_dir = step.dir;

        let r = length(new_pos);
        if (r < r_plus) {
            break;
        }
        if (r > escape_r) {
            let world_dir = normalize(rot_x(new_dir, uniforms.disk_tilt));
            var bg = vec3<f32>(0.0);
            bg += star_color(world_dir, uniforms.star_intensity);
            if (uniforms.skybox_intensity > 0.0) {
                bg += skybox_color(world_dir) * uniforms.skybox_intensity;
            }
            accum_color += (1.0 - accum_alpha) * bg;
            accum_alpha = 1.0;
            break;
        }

        // --- volumetric disk ---
        if (uniforms.disk_quality == 0u) {
            // Off tier: zero-thickness single midplane sample, fixed alpha.
            if (disk_hit(prev, new_pos)) {
                let ty = prev.y / (prev.y - new_pos.y);
                let hit = mix(prev, new_pos, vec3<f32>(ty));
                let s = disk_color_flat(hit, new_dir);
                accum_color += (1.0 - accum_alpha) * s.color * s.density;
                accum_alpha += (1.0 - accum_alpha) * s.density;
                if (accum_alpha > 0.99) { break; }
            }
        } else {
            // Volumetric tier: clip this step to the disk slab and integrate N
            // samples along the clipped segment, weighted by the analytic
            // in-slab length (NOT the RK45 step size). See integrate_disk_segment
            // for why the step-size-decoupled weight is what kills the radial
            // spokes that the previous `* step_len` form produced.
            integrate_disk_segment(prev, new_pos, new_dir, &accum_color, &accum_alpha);
            if (accum_alpha > 0.99) { break; }
        }

        // --- relativistic jets (along the spin axis) ---
        if (uniforms.jets_enabled != 0u) {
            sample_jets(new_pos, new_dir, r_plus, length(new_pos - prev), &accum_color, &accum_alpha);
            if (accum_alpha > 0.99) { break; }
        }

        let ph = planet_hit(prev, new_pos, new_dir);
        if (ph.w > 0.0) {
            accum_color += (1.0 - accum_alpha) * ph.xyz * ph.w;
            accum_alpha += (1.0 - accum_alpha) * ph.w;
            if (accum_alpha > 0.99) { break; }
        }

        if (uniforms.grid_enabled != 0u) {
            let g = grid_hit(prev, new_pos);
            if (g.x + g.y + g.z > 0.0) {
                accum_color += g;
            }
        }

        prev = new_pos;
        pos = new_pos;
        d = new_dir;
    }

    return accum_color;
}

@fragment
fn fragment(in: VertexOutput) -> @location(0) vec4<f32> {
    // Seed the per-pixel hash used by the slab sub-sample jitter (Fix B) and
    // the supersample jitter (Fix A). @builtin(position) is the fragment's
    // framebuffer coordinate; using it makes the noise spatially stable.
    pixel_seed = in.position.xy;

    let aspect = uniforms.resolution.x / uniforms.resolution.y;
    let base_uv = vec2<f32>((in.uv * 2.0 - 1.0).x * aspect, (in.uv * 2.0 - 1.0).y);

    // Per-pixel stochastic supersampling (Fix A). The concentric rings that
    // appear on the disk are higher-order gravitationally-lensed disk images:
    // rays near the critical impact parameter wrap around the hole, and each
    // wrap crosses the disk plane once, projecting to a ring. The integrator
    // renders each wrap as a sharp discrete band rather than a smooth image, so
    // the rings read as an artifact. Firing several jittered sub-rays per pixel
    // and averaging antialiases those sharp band edges into a smooth gradient
    // (the physical Einstein-ring structure is preserved; only its staircase
    // discretization is removed).
    //
    // MAX_AA_SAMPLES is a compile-time cap (WGSL needs static loop bounds);
    // the runtime count comes from uniforms.aa_samples with an early break —
    // the same MAX-with-break form fbm3/ridged_fbm use for dynamic octaves.
    const MAX_AA_SAMPLES = 4u;
    var color_sum = vec3<f32>(0.0);
    var n_done = 0u;
    for (var s: u32 = 0u; s < MAX_AA_SAMPLES; s = s + 1u) {
        if (s >= uniforms.aa_samples) { break; }
        var uv = base_uv;
        if (uniforms.aa_samples > 1u) {
            // ~half-pixel jitter in NDC, seeded by pixel + sample index. The
            // seed is pixel-only (no time) so the pattern is temporally stable.
            let j = hash33(vec3<f32>(pixel_seed, f32(s))) - 0.5;
            uv = uv + j.xy / uniforms.resolution;
        }
        color_sum = color_sum + march(ray_direction(uv));
        n_done = n_done + 1u;
    }

    return vec4<f32>(color_sum / f32(n_done), 1.0);
}