sim.c

// Window-agnostic simulation core. See sim.h for the front-end contract.
// This translation unit emits the glad / dropt / inih single-header impls, so it
// must be linked into every binary exactly once (and the front-ends must not
// re-emit them).

#include "sim.h"
#include "sim_def.h"

#include <ctype.h> // toupper, for the per-sim --help section titles
#include <math.h>
#include <signal.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define DROPT_IMPLEMENTATION // header-only dropt: this is the single TU that emits the impl
#include <dropt.h>           // cli option parser (lib/dropt/)

#define INI_IMPLEMENTATION // header-only inih: this is the single TU that emits the impl
#include <ini.h>           // inih: .ini config parser (lib/ini.h)
#include <limits.h>        // INT_MAX/INT_MIN for parse_int range check
#include <string.h>
#include <strings.h> // strcasecmp

#include "default_params.h" // generated by nob.c from default-params: DEFAULT_PARAMS[]

// Declarations only (no NOB_IMPLEMENTATION): this is for nob_da_append on the GL-name pools.
// The nob_ prefix is kept -- NOB_STRIP_PREFIX would drop Cmd/String_Builder/temp_sprintf/shift
// and friends into this file's namespace.
// Small initial capacity (vs nob's 256) so the pools genuinely realloc while claiming a
// handful of slots. Buffers point at the pool structs rather than into their items precisely
// so that is safe -- at 256 the growth path would never run and the guarantee would rot untested.
#define NOB_DA_INIT_CAP 4
#include <nob.h>

#define GLAD_GL_IMPLEMENTATION // header-only glad: this is the single TU that emits the impl
#include <gl.h>

#include <cglm/cglm.h>
#include <cglm/noise.h>

// width/height are the INTERNAL render (sim) resolution -- logical points divided
// by render_scale. Everything (physics-independent visual passes, density/trail
// buffers, window_shape) is computed against this, then nearest-upscaled to the
// actual framebuffer (window_width/height). This decouples fill cost from device
// DPI (retina renders at 1x) and gives a single perf/quality knob via render_scale.
int width = 0;
int height = 0;
int window_width = 800; // actual framebuffer size in device px; the upscale target
int window_height = 600;
// upscale source-crop (uv), pushed by sim_present. Identity = straight upscale, the
// default for every front-end; macwp -m=-2 narrows it to a centred sub-rect so a
// shared field presents 1:1 on each monitor (crop, not stretch) -- see main-macwp.m.
float present_uv_min[2] = {0.0f, 0.0f};
float present_uv_max[2] = {1.0f, 1.0f};
double render_scale = 0; // internal res = logical points / render_scale (>1 = lower res, faster); --core.render-scale
// SIM_SCALE fixes the simulation coordinate space at logical/SIM_SCALE px, INDEPENDENT of render_scale.
// This makes render_scale a pure fidelity knob (changing -r never alters the dynamics) and a window
// resize just grows/shrinks the field. Set to the historical default render_scale so the default look
// is byte-for-byte preserved; all sim length constants are interpreted in this space.
const double SIM_SCALE = 2.0;
int sim_width = 0, sim_height = 0; // simulation field size in sim px (logical/SIM_SCALE), set in sim_set_dims
float mouse_scale = 1.0f;          // logical-point -> sim-space scale
const char *title = "Pixel Goo";
// runtime-tunable via cli (see parse_args). NOTE: the default VALUES live in the
// default-params file (baked in at build, applied first in parse_args), not here.
// flag-backed globals below init to 0 purely as storage -- default-params is the
// single source of truth and overwrites them before the sim starts. (P is the
// exception: it has no default-params entry, so its fallback constant stays.)
bool fullscreen = true;
bool vsync = true;
bool profile = false;               // --window.profile: glFinish per pass + print ms (disables pipelining)
int max_iterations = 0;             // -N: exit after N PRESENTED frames (0 = unlimited); for benchmarking
int warmup = 0;                     // --core.warmup: simulate this many frames before presenting (skip the boring ramp)
int fps_cap = 0;                    // --window.fps: throttle the average fps to this (0 = uncapped)
bool no_keyfocus_steal = false;     // --window.no-keyfocus-steal: show window w/out grabbing key focus
bool no_border = false;             // --window.no-border: hide the title bar / border in windowed mode
bool no_mouse = false;              // --core.no-mouse: disable the mouse repel (park the cursor far off)
bool mouse_debug = false;           // --core.mouse-debug: draw green dot + trail at cursor
bool corners_debug = false;         // --window.corners-debug: draw green squares in the window corners
bool exclusions_debug = false;      // --core.exclusions-debug: tint excluded regions green
bool edge_debug = false;            // --core.edge-debug: tint the bounding-edge repel bands green
char *dump_path = NULL;             // --window.dump <prefix>: debug. at exit, write frame + density + trail + repel to <prefix>_*.ppm
const char *dump_params_fmt = NULL; // --dump-params[=FORMAT]: print resolved params and exit (NULL = don't)
char *headless_path = NULL;         // --window.headless <out>: pipe raw render-res frames to ffmpeg, no window present
FILE *ffmpeg_pipe = NULL;           // popen handle for the headless encode
unsigned char *headless_px = NULL;  // reusable readback buffer (width*height*3)
int rng_seed = 0;                   // --core.seed <N>: fixed RNG seed for reproducible runs (0 = time-based)
double init_warp = 0;               // --core.init-warp: magnitude of the perlin perturbation on initial positions (0 = pure uniform)
double init_density = 0;            // --core.density/-d: particles per logical-pixel^2. DEFAULT oracle for P unless -p given.
bool p_given = false;               // whether -p/--core.particles was passed -- if so it wins over --core.density
int whichMonitor = 0;

// Texture + framebuffer pools, shared by every sim (see pool_claim). Slots are handed out at
// setup in claim order, so these indices are only valid once buffer_setup has run.
GLuintArray pool_textures = {0};
GLuintArray pool_framebuffers = {0};
PBindex positionBufferIndex1 = -1;
PBindex positionBufferIndex2 = -1;
PBindex velocityBufferIndex1 = -1;
PBindex velocityBufferIndex2 = -1;
PBindex trailBufferIndex1 = -1;
PBindex trailBufferIndex2 = -1;
PBindex renderBufferIndex = -1; // internal colour target, upscaled to the window
PBindex repelBufferIndex = -1;  // force-potential (r) + dispersion (g); single-buffered, stateless

// Physics buffer dims (square-ish texture holding P particles), set in buffer_setup.
int PB_width = 0;
int PB_height = 0;
// PBOs holding position/velocity: filled each frame by glReadPixels(FBO->PBO) and bound
// as vertex attributes for the point passes (faster than vertex texture fetch on this gpu).
//
// Position is a RING. pbo_pos is filled + bound as the point attribute (loc 0) each frame,
// but it is ALSO CPU-mapped to build the tile-coherent draw order. Mapping it for read is the
// stall: on macos GL the map does a conservative pipeline sync -- it blocks even for a slot
// whose readback landed frames ago, because the buffer is part of vertex-array state -- a
// periodic ~4ms (up to ~15ms) hit every sort-every frames. So: (1) ring the buffer, filling +
// attribute-binding a rotating slot each frame; (2) fence each slot when its readback is issued;
// (3) the sort maps the OLDEST slot with GL_MAP_UNSYNCHRONIZED_BIT (skips the driver sync), but
// only once that slot's fence says the readback has actually landed -- otherwise it skips the
// rebuild and keeps the prior order. The sort is a coherence-only optimisation (any positions
// yield a valid permutation), so a skipped or slightly-stale rebuild never affects correctness.
// Velocity is single-buffered: it is never mapped, only streamed as loc 1.
#define POS_PBO_N 4
GLuint pbo_pos[POS_PBO_N] = {0};
int pos_pbo_ring = 0;              // next ring slot to write (this frame's positions)
GLsync pos_fence[POS_PBO_N] = {0}; // per-slot fence: signalled when that slot's readback lands
GLuint pbo_vel = 0;

Shader screenShader = {.name = "screenShader"};
Shader densityShader = {.name = "densityShader"};
Shader positionShader = {.name = "positionShader"};
Shader velocityShader = {.name = "velocityShader"};
Shader copyShader = {.name = "copyShader"};
Shader trailShader = {.name = "trailShader"};
Shader upscaleShader = {.name = "upscaleShader"};
Shader debugShader = {.name = "debugShader"}; // unified debug overlay: exclusions/edge/mouse (fullscreen fill + markers)
Shader repelShader = {.name = "repelShader"};

// Include shader source files. The per-sim ones live in that sim's TU -- these headers
// define the source strings at file scope, so including one in two TUs is a duplicate symbol.
#include "core/debug.h"
#include "core/density.h"
#include "core/position.h"
#include "core/repel.h"
#include "core/upscale.h"

Buffer trailBuffer = {.name = "Trail Buffer", .textures = &pool_textures, .framebuffers = &pool_framebuffers, .current = 0, .other = 1};
Buffer positionBuffer = {.name = "Position buffer", .textures = &pool_textures, .framebuffers = &pool_framebuffers, .current = 0, .other = 1};
Buffer velocityBuffer = {.name = "Velocity buffer", .textures = &pool_textures, .framebuffers = &pool_framebuffers, .current = 0, .other = 1};
Buffer screenBuffer = {.name = "Screen buffer", .textures = NULL, .framebuffers = NULL, .current = 0, .other = 1};
// Internal render target: the visual passes draw here at render resolution, then
// upscaleShader nearest-samples it onto screenBuffer (the actual window).
Buffer renderBuffer = {.name = "Render buffer", .textures = &pool_textures, .framebuffers = &pool_framebuffers, .current = 0, .other = 1};
// Force-potential (r) + dispersion (g) field: single-buffered and stateless -- fully cleared
// and re-splatted (edges + exclusion rects + mouse) from scratch every sim_step, no ping-pong.
Buffer repelBuffer = {.name = "Repel buffer", .textures = &pool_textures, .framebuffers = &pool_framebuffers, .current = 0, .other = 0};

// Alpha blending of each of the fragments
double densityAlpha = 0; // --goo.dens-alpha: per-particle density deposit
double kernelRadius = 0; // --goo.dens-kernel: density splat radius (sim px)
// Force-field tuning (velocity.frag uniforms). All in sim px / sim units, render_scale-independent.
double densityForce = 0;    // --goo.dens-force: density-gradient repel strength
double trailForce = 0;      // --goo.trail-force: trail-gradient attract strength
double repel = 0;           // --core.repel: overall interaction repel strength (walls + exclusions + mouse)
double dragCoefficient = 0; // --core.drag: quadratic velocity damping, both sims
double mouseDrag = 0;       // --core.mouse-drag: gain of the cursor motion sweep (repel field's ba)
int switchFrames = 120;     // --core.switch-frames: length of a sim transition, in FRAMES
double switchEvery = 0;     // --core.switch-every: auto-cycle cadence in seconds (0 = off)
double reach = 0;           // --core.reach: density/interaction VDI sampling radius (sim px)
double trailReach = 0;      // --goo.trail-reach: trail sampling radius (sim px)
double antiStick = 0;       // --core.anti-stick: tiny cardinal drift unsticking motionless particles in repel zones

// The density/trail buffers are heavily downsampled, lerped physics fields -- not
// the rendered image. Scattering every Nth particle into them (with the per-point
// alpha scaled up to compensate) cuts the overdraw/ROP cost without a visible
// change, since the final screen pass still draws all P points. 1 = no subsample.
int densitySubsample = 0; // --goo.dens-sub
int trailSubsample = 0;   // --goo.trail-sub
int dens_every = 0;       // --goo.dens-every: rebuild density field every N frames (reuse between)
int trail_every = 0;      // --goo.trail-every: deposit a rotating 1/N trail subset per frame (decay every frame)

// Screen render: density-weighted cull (--core.cull), 0 = off .. 1 = maximum thinning.
double cullAmount = 0; // --core.cull

// Screen colormap auto-exposure. See parse_args help for the knobs.
double renderHeadroom = 0.0; // --core.headroom: fixed value (0 = auto-track)
double headroomMargin = 0;   // --core.headroom-margin: tracked-max multiplier
double headroomAttack = 0;   // --core.headroom-attack: EMA alpha when density RISES
double headroomRelease = 0;  // --core.headroom-release: EMA alpha when density FALLS
double renderGamma = 0;      // --core.gamma: colormap curve (1 = linear/punchy, <1 lifts faint)
double alphaSpeed = 0;       // --core.alpha-speed: particle speed that reaches full opacity
bool densityNearest = false; // --goo.dens-nearest: GL_NEAREST density (chunky pixels) vs default GL_LINEAR
double headroomPct = 0;      // --core.headroom-pct: track this density percentile, not the max

// Buffer resolution divisor: the density field is render-res / this.
int densityBufferDownsampling = 0;
int repelBufferDownsampling = 0; // --core.repel-downsample: interaction (repel) buffer resolution divisor
int repel_width = 0;             // computed in buffer_setup, like density_width
int repel_height = 0;

double ditherCoefficient = 0; // --core.dither: density-scaled random kick magnitude, both sims
double ditherDensityGain = 0; // --core.dither-gain: gain on the density-scaled dither
double ditherOrtho = 0;       // --core.dither-ortho: density^2-scaled jitter perpendicular to velocity
bool legacyWedge = false;     // --goo.legacy-wedge: old acos trail-integral heading (drops sign of vy)

// Alpha blending of each of the fragments
double trailIntensity = 0; // --goo.trail-intensity: per-particle trail deposit
double trailTau = 0;       // --goo.trail-tau: trail decay time constant in FRAMES (decays every frame)

double trailRadius = 0;          // --goo.trail-kernel: trail splat radius (sim px)
int trailBufferDownsampling = 0; // trail field resolution = render-res / this. --goo.trail-downsample.
double trailVelocityFloor = 0;   // --goo.trail-floor: min velocity for a particle to deposit trail
int trail_width = 0;
int trail_height = 0;

int P = 200000; // --core.particles

// --core.exclusions "rect(x,y,w,h);...": exclusion primitives in logical px. Count 0 = no exclusion
// (today's behaviour). Core, so every front-end exposes it -- see the note in sim.h.
float exclusionRects[MAX_EXCLUSION_RECTS][4] = {{0}};
int exclusionRectCount = 0;

// per-pass GPU timing (only populated under --window.profile). Core's passes first, then the
// active sim's -- pass_name is filled in by sim_setup once that sim is known.
static const char *const core_pass_names[] = {"rdbk", "sort", "draw", "up"};
const int sim_pass_base = (int)(sizeof core_pass_names / sizeof core_pass_names[0]);
double pass_ms[MAX_PASSES] = {0};
const char *pass_name[MAX_PASSES] = {0};
int pass_count = 0;
double pass_t0 = 0; // running lap mark for LAP(); chains across the SimDef calls

// Every registered simulation. All of their cli flags are merged at startup regardless of which one
// runs, since selecting a sim is a runtime decision. Only goo is unconditional; the others need
// their define (see nob.c) -- so this registry, and nothing else, is what makes a sim reachable.
static const SimDef *const sims[] = {
    &SIM_GOO,
#ifdef GOO_SIM_POLLOCK
    &SIM_POLLOCK,
#endif
#ifdef GOO_SIM_BLANK
    &SIM_BLANK,
#endif
};
#define NUM_SIMS (sizeof sims / sizeof sims[0])

// --help section titles, one per sim. dropt keeps the pointer for as long as the context lives, so
// these outlive parse_args rather than being built on its stack.
static char sim_help_headers[MAX_SIMS][32];

// --core.sim's help text, with the registered names listed. Same lifetime reason as the headers
// above, and filled by the same pass over sims[] in parse_args.
static char sim_names_help[128];

// The --window.profile pass set differs per sim (goo has a trail pass, pollock does not), so the
// labels are assembled from whichever sim is active -- at setup AND on every switch, or the
// columns keep the previous sim's names over the new sim's numbers.
static void sim_build_pass_table(void) {
    pass_count = sim_pass_base + active_sim->n_passes;
    if (pass_count > MAX_PASSES) {
        fprintf(stderr, "goo: sim '%s' declares too many profile passes\n", active_sim->name);
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < sim_pass_base; i++)
        pass_name[i] = core_pass_names[i];
    for (int i = 0; i < active_sim->n_passes; i++)
        pass_name[sim_pass_base + i] = active_sim->pass_names[i];
    for (int i = 0; i < MAX_PASSES; i++)
        pass_ms[i] = 0.0; // a switch must not inherit the previous sim's timings
}

// ---- live sim switching ----
// A switch is not a cut: over --core.switch-frames the physics DUTY-CYCLES between the two sims,
// each frame running wholly one or wholly the other, with the incoming share ramping 0 -> 1.
// Force is integrated twice before it becomes visible position, so the particle itself low-passes
// the alternation into a genuine blend -- no extra buffers, no extra passes, one branch here.
// The schedule is error diffusion rather than a random draw: at a 50% share a random draw gives
// runs of four or five identical frames, which is the stutter the whole idea exists to avoid.
//
// Presentation does NOT duty-cycle. The physics integrates the alternation; the display does not,
// so alternating whole frames between two colormaps would strobe. Both sims draw every transition
// frame instead, weighted by the same ramp.
static const SimDef *transition_from = NULL;
static int transition_frame = 0;    // 0 .. switchFrames
static double transition_acc = 0.0; // error-diffusion accumulator
static volatile sig_atomic_t switch_requested = 0;

bool sim_transitioning(void) { return transition_from != NULL; }

void sim_request_switch(void) { switch_requested = 1; }

static void sim_handle_sigusr1(int sig) {
    (void)sig;
    switch_requested = 1; // only a flag: sim_step picks it up
}

void sim_install_signal_handlers(void) {
    signal(SIGUSR1, sim_handle_sigusr1);
}

// Begin a move to the next registered sim, in registry order. Dropped if one is already running:
// queueing them would let a held-down key stack transitions faster than they can play out.
static void sim_begin_switch(void) {
    if (transition_from)
        return;
    size_t at = 0;
    for (size_t i = 0; i < NUM_SIMS; i++)
        if (sims[i] == active_sim)
            at = i;
    const SimDef *next = sims[(at + 1) % NUM_SIMS];
    if (next == active_sim)
        return; // only one sim registered
    transition_from = active_sim;
    active_sim = next;
    transition_frame = 0;
    transition_acc = 0.0;
    // Per-pair setup, before the incoming sim's first duty frame. The outgoing sim's fields are
    // still live and untouched at this point, which is what makes deriving from them possible.
    if (active_sim->enter_from)
        active_sim->enter_from(transition_from);
    sim_build_pass_table(); // the incoming sim's pass set, not the outgoing one's
    fprintf(stdout, "sim: %s -> %s over %d frames\n", transition_from->name, active_sim->name, switchFrames);
}

// The simulation currently running. --core.sim picks the starting one; tab, SIGUSR1 and
// --core.switch-every move it along the registry above.
const SimDef *active_sim = &SIM_GOO;

double gl_refresh_seconds = 30 * 60; // --core.gl-refresh: GL-context recreate cadence (default 30m; 0 = off); wallpaper only

// When set (via sim_restore), the next buffer_setup uploads this saved state
// instead of generating a fresh random field. Consumed (cleared) by buffer_setup.
static const SimSnapshot *g_restore = NULL;

// Monotonic seconds; replaces glfwGetTime
double get_time(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9;
}

// Sleep for ms milliseconds (frame pacing). <=0 returns immediately.
void sleep_ms(double ms) {
    if (ms <= 0.0)
        return;
    struct timespec req = {
        .tv_sec = (time_t)(ms / 1000.0),
        .tv_nsec = (long)(fmod(ms, 1000.0) * 1e6),
    };
    nanosleep(&req, NULL);
}

// Uniform random float in [0, 1]; replaces glm::linearRand
static float frand01(void) {
    return (float)rand() / (float)RAND_MAX;
}

// Exponential moving average; seeds on first sample (prev == 0).
double ema(double prev, double sample, double a) {
    return prev == 0.0 ? sample : a * sample + (1.0 - a) * prev;
}

// Debug dump tooling (--window.dump). Read a framebuffer back and write it as a PPM so it can be
// eyeballed as an image. glReadPixels is bottom-up, so flip to top-down on write.
void dump_ppm_rgb(const char *path, int w, int h, GLuint fbo) {
    glBindFramebuffer(GL_READ_FRAMEBUFFER, fbo);
    unsigned char *px = (unsigned char *)malloc((size_t)w * h * 3);
    // tightly-packed rgb24: default GL_PACK_ALIGNMENT=4 pads rows when w*3 % 4 != 0
    // (e.g. w=125 -> 375), shearing the readback. force byte packing.
    glPixelStorei(GL_PACK_ALIGNMENT, 1);
    glReadPixels(0, 0, w, h, GL_RGB, GL_UNSIGNED_BYTE, px);
    FILE *f = fopen(path, "wb");
    if (f) {
        fprintf(f, "P6\n%d %d\n255\n", w, h);
        for (int y = h - 1; y >= 0; y--)
            fwrite(px + (size_t)y * w * 3, 1, (size_t)w * 3, f);
        fclose(f);
        fprintf(stdout, "dumped %s\n", path);
    }
    free(px);
}

// Single-channel R32F field -> grayscale PPM, normalised by its own max (printed too, so
// the actual value range is visible -- important once density goes unbounded/additive).
void dump_ppm_rgb_scaled(const char *path, int w, int h, GLuint fbo) {
    glBindFramebuffer(GL_READ_FRAMEBUFFER, fbo);
    size_t n = (size_t)w * h * 3;
    float *v = (float *)malloc(n * sizeof(float));
    glReadPixels(0, 0, w, h, GL_RGB, GL_FLOAT, v);
    float mx = 0.0f;
    for (size_t i = 0; i < n; i++)
        if (v[i] > mx)
            mx = v[i];
    // One shared divisor across all channels, not per-channel: normalising each independently
    // would erase the relative abundance between them, which is the point of the dump.
    float norm = mx > 0.0f ? mx : 1.0f;
    unsigned char *px = (unsigned char *)malloc(n);
    for (size_t i = 0; i < n; i++) {
        float t = v[i] / norm;
        px[i] = (unsigned char)(255.0f * (t < 0.0f ? 0.0f : (t > 1.0f ? 1.0f : t)));
    }
    FILE *f = fopen(path, "wb");
    if (f) {
        fprintf(f, "P6\n%d %d\n255\n", w, h);
        for (int y = h - 1; y >= 0; y--)
            fwrite(px + (size_t)y * w * 3, 1, (size_t)w * 3, f);
        fclose(f);
        fprintf(stdout, "dumped %s (max value %.4f)\n", path, mx);
    }
    free(px);
    free(v);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
}

void dump_ppm_scalar(const char *path, int w, int h, GLuint fbo) {
    glBindFramebuffer(GL_READ_FRAMEBUFFER, fbo);
    float *v = (float *)malloc((size_t)w * h * sizeof(float));
    glReadPixels(0, 0, w, h, GL_RED, GL_FLOAT, v);
    float mx = 0.0f;
    for (size_t i = 0; i < (size_t)w * h; i++)
        if (v[i] > mx)
            mx = v[i];
    float norm = mx > 0.0f ? mx : 1.0f;
    unsigned char *px = (unsigned char *)malloc((size_t)w * h);
    for (size_t i = 0; i < (size_t)w * h; i++) {
        float t = v[i] / norm;
        px[i] = (unsigned char)(255.0f * (t < 0.0f ? 0.0f : (t > 1.0f ? 1.0f : t)));
    }
    FILE *f = fopen(path, "wb");
    if (f) {
        fprintf(f, "P5\n%d %d\n255\n", w, h);
        for (int y = h - 1; y >= 0; y--)
            fwrite(px + (size_t)y * w, 1, (size_t)w, f);
        fclose(f);
        fprintf(stdout, "dumped %s (max value %.4f)\n", path, mx);
    }
    free(px);
    free(v);
}

// Tile-coherent draw order. The screen pass binds the goo's scattered points into
// random TBDR tiles, which is the bottleneck; drawing them in screen-tile order via
// an index buffer fixes it. The order is rebuilt by an O(P) counting sort over the
// positions read back into pbo_pos (mapped on the CPU). sort_every controls cadence.
GLuint ibo = 0;
unsigned int *sort_idx = NULL; // P entries: particle indices in tile order
int *sort_key = NULL;          // P entries: per-particle tile bucket (scratch)
int *sort_counts = NULL;       // (nbuckets+1) histogram / prefix-sum scratch
const int sort_tile = 32;      // screen tile size in px (~AGX tile granularity)
int sort_every = 0;            // rebuild the order every N frames (1 = every frame)

// Counting-sort the P particle indices by screen tile, using positions from the OLDEST pbo_pos
// ring slot -- mapped GL_MAP_UNSYNCHRONIZED_BIT (no driver sync) once its fence proves the
// readback landed (see the ring note above). Uploads the result into ibo. O(P), two linear
// passes. Out-of-bounds positions are clamped (only affects coherence, never correctness).
static void rebuild_sort_order(void) {
    int oldest = (pos_pbo_ring + 1) % POS_PBO_N; // least-recently written this ring cycle

    // Only read a slot whose readback has actually landed. Poll the fence (timeout 0, no flush --
    // prior frames' swaps already flushed it). Not signalled, or never filled (startup priming),
    // -> skip this rebuild and keep the prior order (the sort is coherence-only).
    GLsync f = pos_fence[oldest];
    if (f == NULL)
        return;
    GLenum st = glClientWaitSync(f, 0, 0);
    if (st != GL_ALREADY_SIGNALED && st != GL_CONDITION_SATISFIED)
        return;

    int tilesW = (width + sort_tile - 1) / sort_tile;
    int tilesH = (height + sort_tile - 1) / sort_tile;
    int nb = tilesW * tilesH;
    float inv_tile = 1.0f / (float)sort_tile;

    size_t map_bytes = (size_t)PB_width * PB_height * 2 * sizeof(float);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, pbo_pos[oldest]);
    float *pos = (float *)glMapBufferRange(GL_PIXEL_PACK_BUFFER, 0, map_bytes, GL_MAP_READ_BIT | GL_MAP_UNSYNCHRONIZED_BIT);
    if (!pos) {
        glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
        return;
    }

    for (int i = 0; i <= nb; i++)
        sort_counts[i] = 0;
    for (int i = 0; i < P; i++) {
        int tx = (int)(pos[2 * i] * inv_tile);
        int ty = (int)(pos[2 * i + 1] * inv_tile);
        tx = tx < 0 ? 0 : (tx >= tilesW ? tilesW - 1 : tx);
        ty = ty < 0 ? 0 : (ty >= tilesH ? tilesH - 1 : ty);
        int b = ty * tilesW + tx;
        sort_key[i] = b;
        sort_counts[b + 1]++;
    }
    for (int i = 0; i < nb; i++)
        sort_counts[i + 1] += sort_counts[i];
    for (int i = 0; i < P; i++)
        sort_idx[sort_counts[sort_key[i]]++] = (unsigned int)i;

    glUnmapBuffer(GL_PIXEL_PACK_BUFFER);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);

    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, (size_t)P * sizeof(unsigned int), sort_idx, GL_DYNAMIC_DRAW);
}

// dropt int handler that accepts scientific notation ("500", "5e2", "1e6", "2e-1").
// dropt's own dropt_handle_int uses strtol and stops at the 'e'; this uses strtod so
// every int flag takes sci notation too. fractional values round to nearest.
dropt_error parse_int(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL || arg[0] == '\0')
        return dropt_error_insufficient_arguments;
    char *end;
    double v = strtod(arg, &end);
    if (end == arg || *end != '\0')
        return dropt_error_mismatch;
    // round to nearest, then range-check: sci notation makes it trivial to write a
    // value outside int range (e.g. 1e20), and casting an out-of-range double to int
    // is UB. fail loudly instead.
    double r = v >= 0 ? floor(v + 0.5) : ceil(v - 0.5);
    if (r > (double)INT_MAX)
        return dropt_error_overflow;
    if (r < (double)INT_MIN)
        return dropt_error_underflow;
    *(int *)dest = (int)r;
    return dropt_error_none;
}

// dropt handler for --core.render-scale: optional val; bare -r enables low-res mode at 2.0.
static dropt_error parse_render_scale(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL) {
        *(double *)dest = 2.0;
        return dropt_error_none;
    }
    return dropt_handle_double(ctx, opt, arg, dest);
}

// dropt handler for --core.cull: optional val; bare --core.cull enables culling at 0.8.
static dropt_error parse_cull(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL) {
        *(double *)dest = 0.8;
        return dropt_error_none;
    }
    return dropt_handle_double(ctx, opt, arg, dest);
}

// dropt handler for --core.seed: optional val; bare --core.seed picks a random seed and prints it.
static dropt_error parse_seed(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL) {
        int s = (int)(time(NULL) & 0x7fffffff);
        fprintf(stdout, "goo: seed %d\n", s);
        *(int *)dest = s;
        return dropt_error_none;
    }
    return parse_int(ctx, opt, arg, dest);
}

// dropt handler for the particle count: same sci-notation parse as parse_int, but
// also flags that -p was given so it wins over --core.density.
static dropt_error parse_count(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    dropt_error e = parse_int(ctx, opt, arg, dest);
    if (e == dropt_error_none)
        p_given = true;
    return e;
}

// --core.gl-refresh: a duration with an optional unit suffix s/m/h/d (bare = seconds).
// e.g. "10s", "30m", "1h", "1d", "0". Stored as seconds (double).
static dropt_error parse_duration(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL || arg[0] == '\0')
        return dropt_error_insufficient_arguments;
    char *end;
    double v = strtod(arg, &end);
    if (end == arg || v < 0)
        return dropt_error_mismatch;
    double mult = 1.0;
    if (*end != '\0') {
        switch (*end) {
        case 's':
            mult = 1.0;
            break;
        case 'm':
            mult = 60.0;
            break;
        case 'h':
            mult = 3600.0;
            break;
        case 'd':
            mult = 86400.0;
            break;
        default:
            return dropt_error_mismatch;
        }
        if (end[1] != '\0') // only a single-char unit suffix is allowed
            return dropt_error_mismatch;
    }
    *(double *)dest = v * mult;
    return dropt_error_none;
}

// --core.density/-d: plain double, but flag that it was explicitly set (for -p/-d exclusion).
static dropt_error parse_density(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL || arg[0] == '\0')
        return dropt_error_insufficient_arguments;
    char *end;
    double v = strtod(arg, &end);
    if (end == arg || *end != '\0')
        return dropt_error_mismatch;
    *(double *)dest = v;
    return dropt_error_none;
}

// --core.exclusions "rect(x,y,w,h);rect(x,y,w,h);...": exclusion primitives to cut from the sim
// domain, in logical px. Primitive-call syntax: a bare name token followed by a parenthesised,
// comma-separated argument list, ';'-separated calls. Hand-rolled split (no strtok -- same
// manual pointer-walk style as token_flag_name below). Only "rect" (4 numeric args: x,y,w,h)
// is implemented today; any other primitive name, or a malformed call (missing/extra parens,
// wrong arg count, non-numeric field, trailing junk), is a hard reject -- same convention as
// any other flag -- not a silent skip. A primitive count over MAX_EXCLUSION_RECTS is clamped
// with a warning instead, since truncating the tail is a much less surprising failure mode
// than refusing to start.
//
// Storage stays a plain float[4] per rect (no type tag/union) since only one primitive kind
// exists, but the dispatch is keyed off the primitive NAME (not off a fixed field count/
// position), so a second primitive kind later is a new `else if` branch, not a rewrite of the
// call-splitting logic above it.
static dropt_error parse_exclusions(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    if (arg == NULL || arg[0] == '\0')
        return dropt_error_insufficient_arguments;

    int count = 0;
    const char *call_start = arg;
    while (*call_start != '\0') {
        const char *call_end = strchr(call_start, ';');
        size_t call_len = call_end ? (size_t)(call_end - call_start) : strlen(call_start);
        if (call_len > 0) { // tolerate a trailing/empty ';' segment
            const char *paren = memchr(call_start, '(', call_len);
            if (paren == NULL || call_start[call_len - 1] != ')')
                return dropt_error_mismatch; // no "name(...)" shape at all
            size_t name_len = (size_t)(paren - call_start);
            const char *args_start = paren + 1;
            size_t args_len = call_len - name_len - 2; // minus "name(" and the trailing ")"

            if (name_len == 4 && strncmp(call_start, "rect", 4) == 0) {
                if (count >= MAX_EXCLUSION_RECTS) {
                    fprintf(stderr, "goo: --core.exclusions: more than %d primitives given, ignoring the rest\n", MAX_EXCLUSION_RECTS);
                    break;
                }
                float v[4];
                const char *field_start = args_start;
                for (int n = 0; n < 4; n++) {
                    char *end;
                    v[n] = strtof(field_start, &end);
                    if (end == field_start || (n < 3 && *end != ','))
                        return dropt_error_mismatch;
                    field_start = (n < 3) ? end + 1 : end;
                }
                if (field_start != args_start + args_len) // trailing junk before the closing ')'
                    return dropt_error_mismatch;
                exclusionRects[count][0] = v[0];
                exclusionRects[count][1] = v[1];
                exclusionRects[count][2] = v[2];
                exclusionRects[count][3] = v[3];
                count++;
            } else {
                return dropt_error_mismatch; // unrecognised primitive name
            }
        }
        call_start += call_len;
        if (*call_start == ';')
            call_start++;
    }
    *(int *)dest = count;
    return dropt_error_none;
}

// --dump-params[=FORMAT]: optional value. Bare flag defaults to "ini". Just records the format
// string; the actual dump + exit happens at the end of parse_args (once all params are resolved).
static dropt_error parse_dump_params(dropt_context *ctx, const dropt_option *opt, const char *arg, void *dest) {
    (void)ctx;
    (void)opt;
    *(const char **)dest = (arg == NULL || arg[0] == '\0') ? "ini" : arg;
    return dropt_error_none;
}

// Config-from-stdin support. `goo -` reads an .ini config from stdin and parses
// it through the SAME dropt options as the cli, so every flag works in either
// place with one validation path. Each `key = value` becomes a `--key=value`
// token; the cli tokens are appended AFTER the config tokens so the cli wins
// (dropt takes the last value). A key present in both warns to stderr.
//
// NOTE: override detection is by long-flag name only. setting a key in config
// and its short alias on the cli (e.g. `particles` in config + `-p` on cli)
// still overrides correctly but is NOT warned about. add a short->long map here
// if that silent case starts mattering.

// tokens (config + borrowed cli) live for the whole process: dropt's string
// handlers (--window.dump, --window.headless) keep pointers into them, so we never free.
struct token_list {
    char **items;
    size_t count, cap;
};

static void tl_push(struct token_list *tl, char *tok) {
    if (tl->count == tl->cap) {
        tl->cap = tl->cap ? tl->cap * 2 : 16;
        tl->items = realloc(tl->items, tl->cap * sizeof *tl->items);
        if (!tl->items) {
            fprintf(stderr, "goo: out of memory building config tokens\n");
            exit(EXIT_FAILURE);
        }
    }
    tl->items[tl->count++] = tok;
}

// inih handler: turn `[section] name = value` into a `--section.name=value` token, the flag
// spelling those keys have on the cli. A key outside any section keeps its bare name, which is
// then simply not a valid flag -- dropt reports it rather than silently ignoring it.
// Returns 1 (inih: nonzero = ok).
static int config_handler(void *user, const char *section, const char *name, const char *value) {
    struct token_list *tl = user;
    bool sectioned = section && section[0];
    size_t s = sectioned ? strlen(section) + 1 : 0;
    size_t n = strlen(name), v = strlen(value);
    char *tok = malloc(2 + s + n + 1 + v + 1); // "--" + "section." + name + "=" + value + NUL
    if (!tok) {
        fprintf(stderr, "goo: out of memory parsing config\n");
        exit(EXIT_FAILURE);
    }
    if (sectioned)
        sprintf(tok, "--%s.%s=%s", section, name, value);
    else
        sprintf(tok, "--%s=%s", name, value);
    tl_push(tl, tok);
    return 1;
}

// slurp all of stdin into a malloc'd NUL-terminated buffer.
static char *slurp_stdin(void) {
    size_t cap = 4096, len = 0;
    char *buf = malloc(cap);
    if (!buf) {
        fprintf(stderr, "goo: out of memory reading stdin\n");
        exit(EXIT_FAILURE);
    }
    for (;;) {
        if (len + 1 >= cap) {
            cap *= 2;
            buf = realloc(buf, cap);
            if (!buf) {
                fprintf(stderr, "goo: out of memory reading stdin\n");
                exit(EXIT_FAILURE);
            }
        }
        size_t got = fread(buf + len, 1, cap - 1 - len, stdin);
        len += got;
        if (got == 0)
            break;
    }
    buf[len] = '\0';
    return buf;
}

// the long-flag name a token sets, or NULL if it's not a "--name[=...]" token.
// writes the name (sans "--", sans "=value") into out (size cap); returns out or NULL.
static const char *token_flag_name(const char *tok, char *out, size_t cap) {
    if (tok[0] != '-' || tok[1] != '-')
        return NULL;
    const char *p = tok + 2, *eq = strchr(p, '=');
    size_t len = eq ? (size_t)(eq - p) : strlen(p);
    if (len == 0 || len >= cap)
        return NULL;
    memcpy(out, p, len);
    out[len] = '\0';
    return out;
}

// Flags that only make sense for a real window (geometry, input, file output).
// goo-wlwp (wallpaper) excludes these from its option table AND skips any matching
// baked-in default, so they neither show in --help nor parse on its cli.
// NOTE: exclusions and the *-debug overlays are core, not window -- they work on the wallpapers
// via the shared sim core (exclusion physics + the sim_step debug overlay pass). corners-debug is
// window: it's drawn by glScissor in the RGFW front-end loops, which the wallpapers don't have.
static bool is_window_only(const char *name) {
    return strncmp(name, "window.", 7) == 0;
}

// Dump every resolved value param as `key = value`, matching the default-params.ini format, by
// walking the (merged) option table and formatting each entry by its handler. Doubles/ints/bools
// only -- string params (--window.dump, --window.headless) and the meta flags (help, dump-params) are skipped,
// as is `particles` (P isn't derived until sim_setup; the resolved oracle is `density`).
static void dump_params_ini(FILE *out, const dropt_option *opts, size_t n) {
    // Grouped under [section] headers with the prefix stripped, so the output is a config the
    // ini reader accepts verbatim: `goo --dump-params | goo -` has to round-trip.
    const char *emitted[16];
    size_t n_emitted = 0;
    for (size_t pass = 0; pass < n; pass++) {
        const char *want = NULL;
        size_t want_len = 0;
        // pick the next section in first-appearance order
        for (size_t i = 0; i < n && want == NULL; i++) {
            const char *ln = opts[i].long_name;
            const char *dot = ln ? strchr(ln, '.') : NULL;
            if (!dot || opts[i].dest == NULL)
                continue;
            size_t len = (size_t)(dot - ln);
            bool done = false;
            for (size_t k = 0; k < n_emitted; k++)
                if (strlen(emitted[k]) == len && strncmp(emitted[k], ln, len) == 0)
                    done = true;
            if (!done) {
                want = ln;
                want_len = len;
            }
        }
        if (want == NULL)
            break; // every section covered
        if (n_emitted >= sizeof emitted / sizeof emitted[0])
            break;
        char *name = malloc(want_len + 1);
        if (!name) {
            fprintf(stderr, "goo: out of memory dumping params\n");
            exit(EXIT_FAILURE);
        }
        memcpy(name, want, want_len);
        name[want_len] = '\0';
        emitted[n_emitted++] = name;

        fprintf(out, "%s[%s]\n", n_emitted > 1 ? "\n" : "", name);
        for (size_t i = 0; i < n; i++) {
            const dropt_option *o = &opts[i];
            if (o->long_name == NULL || o->dest == NULL)
                continue; // section headers / valueless entries
            const char *dot = strchr(o->long_name, '.');
            if (!dot || (size_t)(dot - o->long_name) != want_len || strncmp(o->long_name, name, want_len) != 0)
                continue;
            const char *key = dot + 1;
            dropt_option_handler_func h = o->handler;
            // parse_duration's dest is plain seconds, and a bare number parses back as seconds,
            // so it dumps as a double like the rest.
            if (h == dropt_handle_double || h == parse_density || h == parse_render_scale || h == parse_cull || h == parse_duration)
                fprintf(out, "%s = %g\n", key, *(const double *)o->dest);
            else if (h == parse_int || h == parse_seed)
                fprintf(out, "%s = %d\n", key, *(const int *)o->dest);
            else if (h == dropt_handle_bool)
                fprintf(out, "%s = %d\n", key, (int)*(const dropt_bool *)o->dest);
            else if (strcmp(key, "sim") == 0)
                fprintf(out, "%s = %s\n", key, active_sim->name); // resolved, so the dump round-trips
            else if (h == parse_exclusions) {
                // Re-emit the primitive-call syntax parse_exclusions accepts, so a config with
                // exclusion rects survives a dump/pipe. Skipped entirely when there are none,
                // rather than emitting an empty value the parser would reject.
                if (exclusionRectCount > 0) {
                    fprintf(out, "%s = ", key);
                    for (int r = 0; r < exclusionRectCount; r++)
                        fprintf(out, "%srect(%g,%g,%g,%g)", r ? ";" : "",
                                exclusionRects[r][0], exclusionRects[r][1],
                                exclusionRects[r][2], exclusionRects[r][3]);
                    fprintf(out, "\n");
                }
            }
            // parse_count (particles), string handlers, etc. -> skipped
        }
    }
    for (size_t k = 0; k < n_emitted; k++)
        free((void *)emitted[k]);
}

// Parse command-line options into the runtime globals above. Exits the process
// on --help or on a bad option; returns normally to let the sim start.
void parse_args(int argc, char **argv, bool wlwp, bool macwp) {
    dropt_bool show_help = 0;
#ifdef GOO_VERSION
    dropt_bool show_version = 0; // --version: only compiled in when GOO_VERSION is baked (CI tag builds)
#endif
    dropt_bool windowed = 0;
    dropt_bool no_vsync = 0;
    dropt_bool prof = 0;
// __APPLE__, not RGFW_MACOS: that's defined by RGFW.h, which this file never includes.
#ifdef __APPLE__
    dropt_bool no_focus = 0;
#endif
    dropt_bool nomouse = 0;
    dropt_bool core_dens_near = 0;
    char *sim_name = NULL;

    // Names for --core.sim's help, from the registry. Each snprintf is guarded on there being room
    // left, since snprintf reports the length it WANTED and that can run past the buffer.
    size_t at_help = 0;
    at_help += (size_t)snprintf(sim_names_help, sizeof sim_names_help, "Which simulation to run (");
    for (size_t i = 0; i < NUM_SIMS && at_help < sizeof sim_names_help; i++)
        at_help += (size_t)snprintf(sim_names_help + at_help, sizeof sim_names_help - at_help,
                                    "%s%s", i ? ", " : "", sims[i]->name);
    if (at_help < sizeof sim_names_help)
        snprintf(sim_names_help + at_help, sizeof sim_names_help - at_help, ").");

    // Core options: the shared simulation, exposed by every front-end. A sim's own flags
    // live in its SimDef and are merged in below; window/wallpaper flags come from the
    // flavour table.
    dropt_option core_opts[] = {
        {'h', "help", "Show this help and exit.", NULL,
         dropt_handle_bool, &show_help, dropt_attr_halt},
#ifdef GOO_VERSION
        {'\0', "version", "Print the version and exit.", NULL,
         dropt_handle_bool, &show_version, dropt_attr_halt},
#endif

        {'\0', NULL, "\nPARTICLES", NULL, NULL, NULL},
        {'p', "core.particles", "Number of particles to simulate (plain or scientific notation, e.g. 200000, 1e6).", "N",
         parse_count, &P},
        {'d', "core.density", "Particles per logical-pixel^2 (sim-area + render-scale independent; default 0.015; excludes -p).", "F",
         parse_density, &init_density},
        {'\0', "core.init-warp", "Magnitude of the perlin warp on initial positions (0 = pure uniform; default 0.13).", "F",
         dropt_handle_double, &init_warp},
        {'\0', "core.repel", "Interaction repel strength: walls + exclusions + mouse (default 0.11).", "F",
         dropt_handle_double, &repel},
        {'\0', "core.drag", "Quadratic velocity damping, both sims (higher = more damping; unrelated to --core.mouse-drag; default 0.19).", "F",
         dropt_handle_double, &dragCoefficient},
        {'\0', "core.mouse-drag", "Gain of the cursor motion sweep (drag).", "F",
         dropt_handle_double, &mouseDrag},
        {'\0', "core.switch-frames", "Length of a sim transition, in frames (physics duty-cycles across it).", "N",
         parse_int, &switchFrames},
        {'\0', "core.switch-every", "Auto-cycle to the next sim every <dur> (s/m/h; 0 = never).", "DUR",
         parse_duration, &switchEvery},
        {'\0', NULL, "\nDENSITY", NULL, NULL, NULL},
        {'\0', "core.dens-force", "Density-gradient force strength (default 0.005).", "F",
         dropt_handle_double, &densityForce},
        {'\0', "core.dens-reach", "Density force sampling radius in sim px (default 20).", "F",
         dropt_handle_double, &reach},
        {'\0', "core.dens-alpha", "Per-particle density deposit (default 0.005).", "F",
         dropt_handle_double, &densityAlpha},
        {'\0', "core.dens-kernel", "Density splat radius in sim px (default 30).", "F",
         dropt_handle_double, &kernelRadius},
        {'\0', "core.dens-sub", "Scatter every Nth particle into EACH density channel (higher = less overdraw).", "N",
         parse_int, &densitySubsample},
        {'\0', "core.dens-downsample", "Density field resolution divisor (render-res / N; higher = coarser/cheaper).", "N",
         parse_int, &densityBufferDownsampling},
        {'\0', "core.dens-every", "Rebuild the density field every N frames (reuse between; cheaper).", "N",
         parse_int, &dens_every},
        {'\0', "core.dens-nearest", "Sample the density field with GL_NEAREST (chunky pixels) instead of GL_LINEAR.", NULL,
         dropt_handle_bool, &core_dens_near},

        // Dither scales with the local density field, so it reads after the DENSITY block above.
        {'\0', NULL, "\nDITHER", NULL, NULL, NULL},
        {'\0', "core.dither", "Density-scaled random kick magnitude, both sims (default 0.025).", "F",
         dropt_handle_double, &ditherCoefficient},
        {'\0', "core.dither-gain", "Gain on the density-scaled dither (>1 = more jitter in dense regions; default 8.5).", "F",
         dropt_handle_double, &ditherDensityGain},
        {'\0', "core.dither-ortho", "Density^2-scaled jitter perpendicular to velocity (default 0.05).", "F",
         dropt_handle_double, &ditherOrtho},

        {'\0', NULL, "\nMISC", NULL, NULL, NULL},
        {'\0', "core.anti-stick", "Anti-stick drift: tiny per-particle cardinal force freeing motionless particles stuck in repel zones (default 7).", "F",
         dropt_handle_double, &antiStick},

        {'\0', NULL, "\nCOLORMAP", NULL, NULL, NULL},
        {'\0', "core.gamma", "Colormap curve: 1 = linear/punchy, <1 lifts faint regions.", "F",
         dropt_handle_double, &renderGamma},
        {'\0', "core.alpha-speed", "Particle speed that reaches full opacity, both sims (default 1.0).", "F",
         dropt_handle_double, &alphaSpeed},
        {'\0', "core.cull", "Screen density cull, 0..1 (bare --core.cull = 0.8; thins denser regions).", "F",
         parse_cull, &cullAmount, dropt_attr_optional_val},
        {'\0', "core.headroom", "Fixed colormap headroom (0 = auto-track the density max).", "F",
         dropt_handle_double, &renderHeadroom},
        {'\0', "core.headroom-margin", "Auto-headroom: multiplier on the tracked max (>1 darker, <1 brighter).", "F",
         dropt_handle_double, &headroomMargin},
        {'\0', "core.headroom-attack", "Auto-headroom: EMA speed when density rises / closes (0..1).", "F",
         dropt_handle_double, &headroomAttack},
        {'\0', "core.headroom-release", "Auto-headroom: EMA speed when density falls / opens (0..1; lower = brightens slower).", "F",
         dropt_handle_double, &headroomRelease},
        {'\0', "core.headroom-pct", "Auto-headroom: track this density percentile (1 = max; <1 ignores hot-spots).", "F",
         dropt_handle_double, &headroomPct},

        {'\0', NULL, "\nDENSITY / TRAIL", NULL, NULL, NULL},
        {'\0', "core.repel-downsample", "Repel (interaction) buffer resolution divisor (render-res / N; higher = coarser edges/exclusions/mouse).", "N",
         parse_int, &repelBufferDownsampling},

        {'\0', NULL, "\nMISC", NULL, NULL, NULL},
        {'r', "core.render-scale", "Internal render resolution divisor (>1 = lower res, faster, chunkier; bare -r = 2).", "N",
         parse_render_scale, &render_scale, dropt_attr_optional_val},
        {'\0', "core.seed", "Fixed RNG seed (bare --core.seed = random but printed for replay; 0 = silent random).", "N",
         parse_seed, &rng_seed, dropt_attr_optional_val},
        {'\0', "core.sort-every", "Rebuild tile-coherent draw order every N frames (0 = never).", "N",
         parse_int, &sort_every},
        {'\0', "core.warmup", "Simulate N frames before presenting -- skips the boring startup ramp (both windowed + headless).", "N",
         parse_int, &warmup},
        {'\0', "dump-params", "Print the resolved params (bare = 'ini' format) and exit.", "FORMAT",
         parse_dump_params, &dump_params_fmt, dropt_attr_optional_val},
        {'\0', "window.no-vsync", "Disable vsync (uncap the frame rate).", NULL,
         dropt_handle_bool, &no_vsync},
        {'\0', "core.no-mouse", "Disable the mouse repel interaction.", NULL,
         dropt_handle_bool, &nomouse},
        {'N', "core.iterations", "Exit after N presented frames (0 = unlimited); for benchmarking.", "N",
         parse_int, &max_iterations},

        // Shared by the window AND the wallpapers -- exclusion physics + the debug overlays live
        // in the sim core (sim_step), so every front-end exposes them.
        // Built from the registry, not spelled out: which sims exist is a compile-time decision
        // (see GOO_SIM_BLANK), so a literal list would go stale on any build but the default one.
        {'\0', "core.sim", sim_names_help, "NAME",
         dropt_handle_string, &sim_name},

        {'\0', NULL, "\nEXCLUSIONS", NULL, NULL, NULL},
        {'\0', "core.exclusions", "Exclusion primitives cut from the sim domain: \"rect(x,y,w,h);rect(x,y,w,h);...\" in logical px (max 16; only 'rect' is implemented).", "PRIMITIVES",
         parse_exclusions, &exclusionRectCount},

    };

    // Core's debug overlays, in their own table so they can be appended AFTER the per-sim blocks:
    // --help prints in table order, and burying the sims behind a DEBUG heading reads as though
    // they were part of it.
    dropt_option core_debug_opts[] = {
        {'\0', NULL, "\nDEBUG", NULL, NULL, NULL},
        {'\0', "core.mouse-debug", "Draw the mouse cursor dot + motion line + repel-radius disc.", NULL,
         dropt_handle_bool, &mouse_debug},
        {'\0', "core.exclusions-debug", "Blend a green tint over excluded regions (debug: visualise exclusion-primitive placement).", NULL,
         dropt_handle_bool, &exclusions_debug},
        {'\0', "core.edge-debug", "Blend a green tint over the bounding-edge repel bands (debug: visualise the edge repel regions).", NULL,
         dropt_handle_bool, &edge_debug},
    };

    // Window-only options: geometry, input, file output -- meaningless for a
    // wallpaper. Appended only for the windowed front-end (see is_window_only).
    dropt_option win_opts[] = {
        {'\0', NULL, "\nWINDOW", NULL, NULL, NULL},
        {'m', "window.monitor", "Index of the monitor to fullscreen onto.", "N",
         parse_int, &whichMonitor},
        {'\0', "window.windowed", "Run in a window instead of fullscreen.", NULL,
         dropt_handle_bool, &windowed},
        {'\0', "window.no-border", "Hide the window border/title bar.", NULL,
         dropt_handle_bool, &no_border},
        {'\0', "window.width", "Window width in pixels (windowed mode).", "N",
         parse_int, &width},
        {'\0', "window.height", "Window height in pixels (windowed mode).", "N",
         parse_int, &height},
        {'\0', "window.fps", "Throttle the average frame rate to N fps (0 = uncapped).", "N",
         parse_int, &fps_cap},

        {'\0', NULL, "\nDEBUG (window-only)", NULL, NULL, NULL},
        {'\0', "window.corners-debug", "Draw green squares in the window corners (debug viewport mapping).", NULL,
         dropt_handle_bool, &corners_debug},
        {'\0', "window.profile", "Print per-pass GPU times in ms (forces glFinish, disables pipelining).", NULL,
         dropt_handle_bool, &prof},
    // NOTE: macos-only -- samespace cover (Spaces, focus-steal) is meaningless on x11/wayland
#ifdef __APPLE__
        {'\0', "window.no-keyfocus-steal", "Show the window without stealing keyboard focus (for benchmarking).", NULL,
         dropt_handle_bool, &no_focus},
#endif
        {'\0', "window.dump", "At exit, write the frame + density, trail and repel fields to <prefix>_*.ppm.", "PREFIX",
         dropt_handle_string, &dump_path},
        {'\0', "window.headless", "Render -N frames at render-res, piped to ffmpeg (e.g. out.mp4); no window. Needs ffmpeg.", "OUT",
         dropt_handle_string, &headless_path},
    };

    // Wallpaper-only options: only meaningful for the wayland layer-shell front-end.
    dropt_option wlwp_opts[] = {
        {'\0', NULL, "\nWALLPAPER", NULL, NULL, NULL},
        {'\0', "core.gl-refresh", "Recreate the GL context every <dur> (s/m/h/d; default 30m; 0 = never) to reclaim the freedreno driver leak. See debug/.", "DUR",
         parse_duration, &gl_refresh_seconds},
    };

    // macos-wallpaper-only options: the wayland wallpaper is single-output, so -m
    // (clone to a monitor / all monitors) only makes sense for goo-macwp.
    dropt_option macwp_opts[] = {
        {'\0', NULL, "\nWALLPAPER", NULL, NULL, NULL},
        {'m', "window.monitor", "Monitor for the wallpaper (-1 = clone+stretch to all; -2 = clone to all, centre-crop no stretch; default 0 = main).", "N",
         parse_int, &whichMonitor},
    };

    // Assemble the final table: core + flavor-specific extras + sentinel.
    // Every registered sim's flags are merged, not just the active one's: a switch is a
    // runtime decision, so all of them have to resolve at startup.
    size_t n_sim = 0;
    for (size_t i = 0; i < NUM_SIMS; i++)
        n_sim += sims[i]->n_opts;
    size_t n_core = sizeof core_opts / sizeof core_opts[0];
    dropt_option *extra = macwp ? macwp_opts : wlwp ? wlwp_opts
                                                    : win_opts;
    size_t n_extra = macwp  ? sizeof macwp_opts / sizeof macwp_opts[0]
                     : wlwp ? sizeof wlwp_opts / sizeof wlwp_opts[0]
                            : sizeof win_opts / sizeof win_opts[0];
    size_t n_core_debug = sizeof core_debug_opts / sizeof core_debug_opts[0];
    // +NUM_SIMS: each sim's block is preceded by its own --help section header. Without one it
    // would print under whichever section the preceding table happened to end on.
    dropt_option *options = malloc((n_core + n_sim + NUM_SIMS + n_core_debug + n_extra + 1) * sizeof(dropt_option));
    if (!options) {
        fprintf(stderr, "goo: out of memory setting up option parser\n");
        exit(EXIT_FAILURE);
    }
    memcpy(options, core_opts, n_core * sizeof(dropt_option));
    size_t at = n_core;
    for (size_t i = 0; i < NUM_SIMS && i < MAX_SIMS; i++) {
        // dropt reads an entry with no long name as a section title; the leading newline spaces
        // it off the previous block, matching the headers in the tables above.
        snprintf(sim_help_headers[i], sizeof sim_help_headers[i], "\n%s", sims[i]->name);
        for (char *c = sim_help_headers[i]; *c; c++)
            *c = (char)toupper((unsigned char)*c);
        options[at++] = (dropt_option){'\0', NULL, sim_help_headers[i], NULL, NULL, NULL};
        memcpy(options + at, sims[i]->opts, sims[i]->n_opts * sizeof(dropt_option));
        at += sims[i]->n_opts;
    }
    memcpy(options + at, core_debug_opts, n_core_debug * sizeof(dropt_option));
    at += n_core_debug;
    memcpy(options + at, extra, n_extra * sizeof(dropt_option));
    memset(options + at + n_extra, 0, sizeof(dropt_option)); // sentinel

    const char *prog = macwp ? "goo-macwp" : wlwp ? "goo-wlwp"
                                                  : "goo";

    dropt_context *ctx = dropt_new_context(options);
    if (ctx == NULL) {
        fprintf(stderr, "goo: out of memory setting up option parser\n");
        exit(EXIT_FAILURE);
    }
    dropt_allow_concatenated_arguments(ctx, true);

    // Build the token stream dropt parses. If a bare `-` appears in argv, read an
    // .ini config from stdin first; its tokens go in front so the cli (appended
    // after) overrides on conflict.
    struct token_list cfg = {0};
    bool want_stdin = false;
    bool want_help = false;
    for (int i = 1; i < argc; i++) {
        if (strcmp(argv[i], "-") == 0)
            want_stdin = true;
        // help wins over everything: honour -h/--help even if other args are bad.
        if (strcmp(argv[i], "-h") == 0 || strcmp(argv[i], "--help") == 0)
            want_help = true;
    }

    // Print help before parsing/stdin so a bad arg -- or a piped `-` we'd block on
    // -- can't suppress it.
    if (want_help) {
        fprintf(stdout, "usage: %s [options]\n", prog);
        dropt_help_params hp;
        dropt_init_help_params(&hp);
        hp.indent = 2;
        hp.description_start_column = 26;
        hp.blank_lines_between_options = false;
        dropt_print_help(stdout, ctx, &hp);
        dropt_free_context(ctx);
        exit(EXIT_SUCCESS);
    }

    // Layer 0: the baked-in defaults from default-params (source of truth, see
    // nob.c). These come first so anything below overrides them silently. For wlwp,
    // skip the window.* defaults -- those flags aren't in a wallpaper's table, so dropt
    // would reject the tokens.
    for (size_t i = 0; i < sizeof DEFAULT_PARAMS / sizeof DEFAULT_PARAMS[0]; i++) {
        char key[128];
        if (wlwp && token_flag_name(DEFAULT_PARAMS[i], key, sizeof key) && is_window_only(key))
            continue;
        tl_push(&cfg, (char *)DEFAULT_PARAMS[i]);
    }

    // Layer 1: stdin config (only when `-` is present).
    size_t stdin_start = cfg.count;
    if (want_stdin) {
        char *buf = slurp_stdin();
        int err = ini_parse_string(buf, config_handler, &cfg);
        free(buf);
        if (err < 0) {
            fprintf(stderr, "goo: cannot read config from stdin\n");
            exit(EXIT_FAILURE);
        }
        if (err > 0) {
            fprintf(stderr, "goo: config parse error on stdin line %d\n", err);
            exit(EXIT_FAILURE);
        }
    }
    size_t stdin_end = cfg.count;

    // Layer 2: the cli, appended last so it wins (dropt keeps the last value).
    for (int i = 1; i < argc; i++) {
        if (strcmp(argv[i], "-") == 0)
            continue; // the stdin sentinel, not a flag
        tl_push(&cfg, argv[i]);
    }

    // Warn for every stdin-config key the cli overrides (long-flag names only;
    // see NOTE). Defaults overridden by stdin/cli are expected and stay silent.
    for (size_t i = stdin_start; i < stdin_end; i++) {
        char ckey[128];
        if (!token_flag_name(cfg.items[i], ckey, sizeof ckey))
            continue;
        for (size_t j = stdin_end; j < cfg.count; j++) {
            char lkey[128];
            if (token_flag_name(cfg.items[j], lkey, sizeof lkey) && strcmp(ckey, lkey) == 0) {
                fprintf(stderr, "goo: config '%s' overridden by cli\n", ckey);
                break;
            }
        }
    }

    dropt_parse(ctx, (int)cfg.count, cfg.items);

    if (dropt_get_error(ctx) != dropt_error_none) {
        fprintf(stderr, "goo: %s\n", dropt_get_error_message(ctx));
        dropt_free_context(ctx);
        exit(EXIT_FAILURE);
    }

#ifdef GOO_VERSION
    if (show_version) {
        fprintf(stdout, "%s\n", GOO_VERSION);
        dropt_free_context(ctx);
        free(options);
        exit(EXIT_SUCCESS);
    }
#endif

    // Resolve --core.sim now that parsing is done. Unknown names are an error rather than a
    // silent fallback: a typo would otherwise run the wrong sim with the right-looking flags.
    if (sim_name) {
        const SimDef *found = NULL;
        for (size_t i = 0; i < NUM_SIMS; i++)
            if (strcmp(sims[i]->name, sim_name) == 0)
                found = sims[i];
        if (!found) {
            fprintf(stderr, "goo: unknown --core.sim '%s' (have:", sim_name);
            for (size_t i = 0; i < NUM_SIMS; i++)
                fprintf(stderr, " %s", sims[i]->name);
            fprintf(stderr, ")\n");
            exit(EXIT_FAILURE);
        }
        active_sim = found;
    }

    // --dump-params[=FORMAT]: every param is resolved now (and we're still before any window/GL
    // setup, since parse_args runs first). Print + exit. Only 'ini' format so far.
    if (dump_params_fmt != NULL) {
        if (strcmp(dump_params_fmt, "ini") != 0) {
            fprintf(stderr, "goo: unknown --dump-params format '%s' (supported: ini)\n", dump_params_fmt);
            dropt_free_context(ctx);
            free(options);
            exit(EXIT_FAILURE);
        }
        dump_params_ini(stdout, options, n_core + n_sim + n_extra);
        dropt_free_context(ctx);
        free(options);
        exit(EXIT_SUCCESS);
    }

    dropt_free_context(ctx);
    free(options);

    // Apply the negating switches onto the defaults.
    if (windowed)
        fullscreen = false;
    if (no_vsync)
        vsync = false;
    if (prof)
        profile = true;
#ifdef __APPLE__
    if (no_focus)
        no_keyfocus_steal = true;
#endif
    for (size_t i = 0; i < NUM_SIMS; i++)
        if (sims[i]->apply_opts)
            sims[i]->apply_opts();
    if (nomouse)
        no_mouse = true;
    if (core_dens_near)
        densityNearest = true;

    // --core.dens-sub is per channel, so a 3-channel field wants 3x the points at the same
    // setting -- below the channel count that would need more points than there are particles,
    // and density_points clamps instead. Say so once rather than silently under-depositing.
    if (densitySubsample > 0 && densitySubsample < 3)
        fprintf(stderr, "goo: --core.dens-sub=%d is below the 3 channels a multi-channel sim uses;"
                        " its field will be under-deposited (clamped at one pass over every particle)\n",
                densitySubsample);

    // Guard against values that would crash or hang the sim.
    if (P < 1) {
        fprintf(stderr, "goo: --core.particles must be >= 1\n");
        exit(EXIT_FAILURE);
    }
    // wlwp gets width/height from the layer surface (set after parse), so they're
    // legitimately 0 here -- only the windowed front-end takes them from flags.
    if (!wlwp && (width < 1 || height < 1)) {
        fprintf(stderr, "goo: --window.width/--window.height must be >= 1\n");
        exit(EXIT_FAILURE);
    }
    if (render_scale < 1.0) {
        fprintf(stderr, "goo: --core.render-scale must be >= 1\n");
        exit(EXIT_FAILURE);
    }
    if (densityBufferDownsampling < 1 || trailBufferDownsampling < 1 || repelBufferDownsampling < 1) {
        fprintf(stderr, "goo: --goo.dens-downsample/--goo.trail-downsample/--core.repel-downsample must be >= 1\n");
        exit(EXIT_FAILURE);
    }
    // Decay time constant feeds exp(-1/tau): tau <= 0 gives alpha 0 (trail wiped every frame),
    // which silently kills the trail dynamics. Reject it.
    if (trailTau <= 0.0) {
        fprintf(stderr, "goo: --goo.trail-tau must be > 0\n");
        exit(EXIT_FAILURE);
    }
    if (densitySubsample < 1 || trailSubsample < 1) {
        fprintf(stderr, "goo: --goo.dens-sub/--goo.trail-sub must be >= 1\n");
        exit(EXIT_FAILURE);
    }
    if (warmup < 0) {
        fprintf(stderr, "goo: --core.warmup must be >= 0\n");
        exit(EXIT_FAILURE);
    }
}

//=====================================================
//
//   ####  ##   ##    ###    ####    #####  #####
//  ##     ##   ##   ## ##   ##  ##  ##     ##  ##
//   ###   #######  ##   ##  ##  ##  #####  #####
//     ##  ##   ##  #######  ##  ##  ##     ##  ##
//  ####   ##   ##  ##   ##  ####    #####  ##   ##
//
//=====================================================

void shader_setup(void) {
    fprintf(stdout, "Compiling shaders...\n");
    // Every registered sim is compiled and allocated, not just the active one: a switch has to
    // be a pointer change, and a transition needs both sims' fields live at once.
    for (size_t i = 0; i < NUM_SIMS; i++)
        sims[i]->shaders();

    // Position integration is shared: every sim advances positions from velocities the same
    // way, so core owns the program and each sim drives it from its own integrate().
    shader_create(&positionShader);
    shader_compile(&positionShader, GL_VERTEX_SHADER, core_position_VertexShaderSource);
    shader_compile(&positionShader, GL_FRAGMENT_SHADER, core_position_FragmentShaderSource);
    shader_link(&positionShader);

    // Density is a core service now: one program, one splat, every sim's field built the same way.
    shader_create(&densityShader);
    shader_compile(&densityShader, GL_VERTEX_SHADER, core_density_VertexShaderSource);
    shader_compile(&densityShader, GL_FRAGMENT_SHADER, core_density_FragmentShaderSource);
    shader_link(&densityShader);

    shader_create(&upscaleShader);
    shader_compile(&upscaleShader, GL_VERTEX_SHADER, core_upscale_VertexShaderSource);
    shader_compile(&upscaleShader, GL_FRAGMENT_SHADER, core_upscale_FragmentShaderSource);
    shader_link(&upscaleShader);

    shader_create(&debugShader);
    shader_compile(&debugShader, GL_VERTEX_SHADER, core_debug_VertexShaderSource);
    shader_compile(&debugShader, GL_FRAGMENT_SHADER, core_debug_FragmentShaderSource);
    shader_link(&debugShader);

    shader_create(&repelShader);
    shader_compile(&repelShader, GL_VERTEX_SHADER, core_repel_VertexShaderSource);
    shader_compile(&repelShader, GL_FRAGMENT_SHADER, core_repel_FragmentShaderSource);
    shader_link(&repelShader);

    // debugShader (created above) must be linked before updateShaderWindowShape ->
    // push_exclusion_rects (below), which pushes the render-space exclusion rects to it.
    updateShaderWindowShape(width, height);
}

// Push the exclusion-rect list (--core.exclusions) to every shader that consumes it, each converted
// into that shader's own coordinate space. Called from updateShaderWindowShape so it's
// refreshed on resize too -- both conversions below depend on the current window size.
static void push_exclusion_rects(void) {
    int count = exclusionRectCount > MAX_EXCLUSION_RECTS ? MAX_EXCLUSION_RECTS : exclusionRectCount;
    shader_set_uniform_int(&repelShader, "exclusion_rect_count", count);
    shader_set_uniform_int(&debugShader, "exclusion_rect_count", count);
    for (int i = 0; i < MAX_EXCLUSION_RECTS; i++) {
        char name[32];
        snprintf(name, sizeof name, "exclusion_rects[%d]", i);
        float x = exclusionRects[i][0], y = exclusionRects[i][1], w = exclusionRects[i][2], h = exclusionRects[i][3];

        // Sim-space bounds (repel.frag's wall splat): logical px -> sim px is a plain
        // scale (mouse_scale) -- sim-space y is top-down, same as logical px, no flip needed
        // (matches how mouse_position is scaled in the front-end).
        float sim_rect[4] = {x * mouse_scale, y * mouse_scale, (x + w) * mouse_scale, (y + h) * mouse_scale};
        shader_set_uniform_vec(&repelShader, name, 4, sim_rect);

        // Render-buffer-pixel bounds (debug.frag's overlay): logical px -> render px is a
        // render_scale divide PLUS a y-flip (gl_FragCoord is y-up/bottom-origin; logical px,
        // like the mouse, is top-down) -- mirrors how screen.vert/density.vert project
        // positions into their own render targets.
        float render_rect[4] = {
            x / (float)render_scale,
            (float)height - (y + h) / (float)render_scale,
            (x + w) / (float)render_scale,
            (float)height - y / (float)render_scale,
        };
        shader_set_uniform_vec(&debugShader, name, 4, render_rect);
    }
}

void updateShaderWindowShape(int new_width, int new_height) {
    (void)new_width;
    (void)new_height; // shape comes from the global sim/render dims now

    // window_shape is the SIM space (logical/SIM_SCALE): all physics, density/trail splat, and the
    // screen vertex pass normalise particle positions by it, so the dynamics are render_scale-free.
    float sim_shape[2] = {(float)sim_width, (float)sim_height};
    shader_set_uniform_vec(&positionShader, "window_shape", 2, sim_shape);
    shader_set_uniform_vec(&repelShader, "window_shape", 2, sim_shape);
    shader_set_uniform_vec(&densityShader, "window_shape", 2, sim_shape);

    // Render-buffer pixels, for the passes that sample by gl_FragCoord rather than sim position.
    float render_shape[2] = {(float)width, (float)height};
    for (size_t i = 0; i < NUM_SIMS; i++)
        sims[i]->update_shapes(sim_shape, render_shape);

    push_exclusion_rects();

    // The screen buffer is the final upscale target, so it takes the actual window (framebuffer) size.
    screenBuffer.width = window_width;
    screenBuffer.height = window_height;
}

void updateShaderBufferShape(int new_width, int new_height) {
    // buffer_size is dead: the point passes now read positions as PBO attributes
    // instead of mapping gl_VertexID -> texel, so no shader references it anymore.
    (void)new_width;
    (void)new_height;
}

//===================================================
//
//  #####   ##   ##  #####  #####  #####  #####
//  ##  ##  ##   ##  ##     ##     ##     ##  ##
//  #####   ##   ##  #####  #####  #####  #####
//  ##  ##  ##   ##  ##     ##     ##     ##  ##
//  #####    #####   ##     ##     #####  ##   ##
//
//===================================================

// ---- auto-exposure follower (see sim.h) ----
// A direct glReadPixels(FBO -> client) would force a full CPU/GPU sync every update, a periodic
// multi-ms stall. Instead the readback goes async into a ring of PBOs and we map the OLDEST slot,
// filled HEADROOM_PBO_N-1 updates ago, so its DMA is long done and the map never stalls. The
// follower is a slow envelope, so a few frames of latency is invisible.

static int cmp_float(const void *a, const void *b) {
    float fa = *(const float *)a, fb = *(const float *)b;
    return (fa > fb) - (fa < fb);
}

// ---- density field ----
// Attribute 2 is the per-particle channel index. It has to be set per PASS, not once at setup:
// every sim is resident at the same time, so a 1-channel sim's splat runs in the same frame as a
// 3-channel sim's and would otherwise inherit whatever the other left bound -- which would send
// two thirds of a single-channel sim's particles into channels its texture does not store.
#define DENSITY_CHANNEL_ATTRIB 2

// Points to draw per rebuild. --core.dens-sub is defined PER CHANNEL: a 3-channel field draws three
// times as many points as a 1-channel one at the same setting, so each channel ends up with the same
// particle count -- and so the same field magnitude -- as a single-channel sim's. Without that, a
// shared --core.dens-force would be dimensionally equal between the sims but not equal in effect.
static int density_points(const Density *d) {
    long want = (long)P * d->channels / (densitySubsample > 0 ? densitySubsample : 1);
    return want > (long)P ? P : (int)want;
}

void density_create(Density *d, const char *name, int channels) {
    d->channels = channels;
    d->width = width / densityBufferDownsampling + 1;
    d->height = height / densityBufferDownsampling + 1;
    d->buffer.name = name;
    d->buffer.textures = &pool_textures;
    d->buffer.framebuffers = &pool_framebuffers;
    d->buffer.other = 0; // single-buffered: the splat clears and rewrites it every rebuild
    d->buffer.minmag_filter = densityNearest ? GL_NEAREST : GL_LINEAR;
    d->buffer.wrap_st = GL_REPEAT;
    d->buffer.dim = channels == 1 ? BE_1D : BE_3D;
    buffer_allocate(&d->buffer, current, pool_claim(1), d->width, d->height, NULL);

    headroom_forget(&d->headroom); // fresh GL context: the ring's old names are meaningless
    headroom_reset(&d->headroom, d->width, d->height, channels);
}

void density_resize(Density *d) {
    d->width = width / densityBufferDownsampling + 1;
    d->height = height / densityBufferDownsampling + 1;
    buffer_reallocate(&d->buffer, current, d->width, d->height);
    headroom_reset(&d->headroom, d->width, d->height, d->channels);
}

void density_splat(Density *d) {
    shader_use(&densityShader);
    shader_set_uniform_int(&densityShader, "density_buffer_downsampling", densityBufferDownsampling);
    shader_set_uniform_float(&densityShader, "density_alpha", densityAlpha * densitySubsample);
    // kernelRadius is in sim px; the density buffer is at render res, so convert sim->render-px.
    shader_set_uniform_float(&densityShader, "kernel_radius", kernelRadius * (float)(SIM_SCALE / render_scale));

    if (d->channels > 1 && d->channel_vbo) {
        glBindBuffer(GL_ARRAY_BUFFER, d->channel_vbo);
        glEnableVertexAttribArray(DENSITY_CHANNEL_ATTRIB);
        glVertexAttribPointer(DENSITY_CHANNEL_ATTRIB, 1, GL_FLOAT, GL_FALSE, 0, 0);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
    } else {
        glDisableVertexAttribArray(DENSITY_CHANNEL_ATTRIB);
        glVertexAttrib1f(DENSITY_CHANNEL_ATTRIB, 0.0f);
    }

    buffer_bind(&d->buffer, current);
    glBlendFunc(GL_ONE, GL_ONE); // additive density accumulation (unbounded)
    buffer_update_n(&d->buffer, density_points(d));
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // restore the global blend

    headroom_update(&d->headroom, &d->buffer);
}

void headroom_forget(HeadroomState *hs) {
    memset(hs->pbo, 0, sizeof hs->pbo);
}

void headroom_reset(HeadroomState *hs, int w, int h, int channels) {
    hs->w = w;
    hs->h = h;
    hs->channels = channels;
    if (hs->pbo[0] == 0)
        glGenBuffers(HEADROOM_PBO_N, hs->pbo); // fresh names per GL context
    size_t bytes = (size_t)w * h * channels * sizeof(float);
    for (int i = 0; i < HEADROOM_PBO_N; i++) {
        glBindBuffer(GL_PIXEL_PACK_BUFFER, hs->pbo[i]);
        glBufferData(GL_PIXEL_PACK_BUFFER, bytes, NULL, GL_STREAM_READ);
    }
    glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
    hs->ring = 0;
    hs->primed = 0; // re-prime: the old ring contents are the wrong size / stale
    // Generously padded so a grow-resize can't overflow it before the next reset.
    hs->readback = (float *)realloc(hs->readback, (size_t)(w + 8) * (h + 8) * channels * sizeof(float));
    // Seed 0 in auto mode so ema()'s seed-on-zero snaps to the FIRST measured level instead of
    // ramping from a guessed constant (that ramp showed as faint -> bright over the first frames).
    hs->ema = (float)(renderHeadroom > 0.0 ? renderHeadroom : 0.0);
}

void headroom_update(HeadroomState *hs, const Buffer *density) {
    if (renderHeadroom > 0.0)
        return; // fixed headroom: nothing to track
    int n = hs->w * hs->h;
    GLenum fmt = hs->channels == 1 ? GL_RED : (hs->channels == 3 ? GL_RGB : GL_RGBA);
    size_t vals = (size_t)n * hs->channels;

    // Async DMA of the fresh density into this ring slot -- glReadPixels into a bound
    // GL_PIXEL_PACK_BUFFER is a GPU-side copy that returns immediately (no CPU sync).
    glBindFramebuffer(GL_READ_FRAMEBUFFER, density->framebuffers->items[density->current]);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, hs->pbo[hs->ring]);
    glReadPixels(0, 0, hs->w, hs->h, fmt, GL_FLOAT, 0);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);

    bool have = false;
    if (hs->primed >= HEADROOM_PBO_N) {
        // Steady state: map the OLDEST slot (its DMA is long done, so the map does not stall).
        int oldest = (hs->ring + 1) % HEADROOM_PBO_N;
        glBindBuffer(GL_PIXEL_PACK_BUFFER, hs->pbo[oldest]);
        void *p = glMapBufferRange(GL_PIXEL_PACK_BUFFER, 0, vals * sizeof(float), GL_MAP_READ_BIT);
        if (p) {
            memcpy(hs->readback, p, vals * sizeof(float));
            glUnmapBuffer(GL_PIXEL_PACK_BUFFER);
            have = true;
        }
        glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
    } else {
        // Startup priming: the ring isn't full yet. One synchronous read so the EMA tracks from
        // the first frames (a few stalls during the warmup ramp), then the async path takes over.
        glBindFramebuffer(GL_READ_FRAMEBUFFER, density->framebuffers->items[density->current]);
        glReadPixels(0, 0, hs->w, hs->h, fmt, GL_FLOAT, hs->readback);
        glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
        have = true;
    }
    hs->ring = (hs->ring + 1) % HEADROOM_PBO_N;
    if (hs->primed < HEADROOM_PBO_N)
        hs->primed++;
    if (!have)
        return;

    // Reduce to one scalar per texel before the percentile -- see the note in sim.h.
    if (hs->channels > 1) {
        for (int i = 0; i < n; i++) {
            float m = hs->readback[(size_t)i * hs->channels];
            for (int c = 1; c < hs->channels; c++) {
                float v = hs->readback[(size_t)i * hs->channels + c];
                if (v > m)
                    m = v;
            }
            hs->readback[i] = m;
        }
    }

    qsort(hs->readback, n, sizeof(float), cmp_float);
    int k = (int)(headroomPct * (n - 1));
    if (k < 0)
        k = 0;
    else if (k >= n)
        k = n - 1;
    float target = hs->readback[k] * (float)headroomMargin;
    // Asymmetric follower: a rising target (scene densifying) closes fast to avoid clipping;
    // a falling one (opening up) brightens slowly so it doesn't pop.
    double a = (hs->ema > 0.0f && target < hs->ema) ? headroomRelease : headroomAttack;
    hs->ema = (float)ema(hs->ema, target, a);
}

// Reserve n consecutive pool slots and generate their GL names. A slot index doubles as the
// texture unit buffer_allocate binds to (glActiveTexture(GL_TEXTURE0 + index)), so the pool is
// bounded by the driver's combined texture-unit limit rather than by an arbitrary constant.
PBindex pool_claim(int n) {
    GLint max_units = 0;
    glGetIntegerv(GL_MAX_COMBINED_TEXTURE_IMAGE_UNITS, &max_units);
    if ((int)pool_textures.count + n > max_units) {
        fprintf(stderr, "goo: out of texture units (%d needed, %d available)\n",
                (int)pool_textures.count + n, max_units);
        exit(EXIT_FAILURE);
    }
    PBindex first = (PBindex)pool_textures.count;
    for (int i = 0; i < n; i++) {
        nob_da_append(&pool_textures, 0);
        nob_da_append(&pool_framebuffers, 0);
    }
    glGenTextures(n, &pool_textures.items[first]);
    glGenFramebuffers(n, &pool_framebuffers.items[first]);
    return first;
}

// Drop the pools without touching GL. buffer_setup uses this: on a --core.gl-refresh the old context
// is already destroyed by the time it runs, so the old names are meaningless and deleting them
// would act on whatever the NEW context happens to have assigned those numbers to.
static void pool_reset(void) {
    nob_da_free(pool_textures);
    nob_da_free(pool_framebuffers);
    pool_textures = (GLuintArray){0};
    pool_framebuffers = (GLuintArray){0};
}

void sim_teardown_gl(void) {
    if (pool_textures.count > 0) {
        glDeleteTextures((GLsizei)pool_textures.count, pool_textures.items);
        glDeleteFramebuffers((GLsizei)pool_framebuffers.count, pool_framebuffers.items);
    }
    pool_reset();
}

void buffer_setup(void) {
    fprintf(stdout, "Setting up buffers...\n");

    // Derived buffer sizes, now that width/height are final (set in the front-end's window setup).
    // The density field's dims live in each sim's Density unit and are computed by density_create.
    repel_width = width / repelBufferDownsampling + 1;
    repel_height = height / repelBufferDownsampling + 1;
    trail_width = width / trailBufferDownsampling + 1;
    trail_height = height / trailBufferDownsampling + 1;

    // Vertex array. The point passes get their attributes from the PBOs (wired below);
    // the fullscreen-quad passes use gl_VertexID -- so the VAO needs no static buffer.
    GLuint vao;
    glGenVertexArrays(1, &vao);
    glBindVertexArray(vao);

    // Calculate size of the physics buffer
    int PBwidth = ceil(sqrt(P));
    int PBheight = ceil(P / sqrt(P));
    PB_width = PBwidth;
    PB_height = PBheight;
    fprintf(stdout, "%d points\n", P);
    // Initial density (inverse of -d): P over logical screen area, so it's comparable across
    // sim sizes / render scales regardless of whether -p or -d set the count.
    double logical_area = (double)width * height * render_scale * render_scale;
    fprintf(stdout, "Initial density: %.4g particles/logical-px^2\n", logical_area > 0 ? P / logical_area : 0.0);
    fprintf(stdout, "Physics framebuffer shape: %d x %d\n", PBwidth, PBheight);
    updateShaderBufferShape(PBwidth, PBheight);

    // PBOs sized to the full physics buffer (RG = 2 floats/texel). Filled each frame
    // by glReadPixels(FBO->PBO), then reused as GL_ARRAY_BUFFER so the point passes
    // read positions/velocities as streamed attributes instead of vertex texture fetch.
    size_t pb_bytes = (size_t)PBwidth * PBheight * 2 * sizeof(float);
    // Position PBO ring (see the ring note by the pbo_pos declaration). P is fixed, so these
    // are never resized -- only regenerated if buffer_setup re-runs on a context recreate.
    glGenBuffers(POS_PBO_N, pbo_pos);
    for (int i = 0; i < POS_PBO_N; i++) {
        glBindBuffer(GL_PIXEL_PACK_BUFFER, pbo_pos[i]);
        glBufferData(GL_PIXEL_PACK_BUFFER, pb_bytes, NULL, GL_DYNAMIC_COPY);
    }
    pos_pbo_ring = 0;
    // Abandon any fences from a prior context (they died with it); a fresh ring starts unfilled.
    for (int i = 0; i < POS_PBO_N; i++)
        pos_fence[i] = 0;
    glGenBuffers(1, &pbo_vel);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, pbo_vel);
    glBufferData(GL_PIXEL_PACK_BUFFER, pb_bytes, NULL, GL_DYNAMIC_COPY);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);

    // Wire the PBOs as vertex attributes (position = loc 0, velocity = loc 1) on the
    // still-bound VAO, so the point passes stream them instead of vertex texture fetch.
    // Position starts on ring slot 0; sim_step re-points it at the freshly filled slot each frame.
    glBindBuffer(GL_ARRAY_BUFFER, pbo_pos[0]);
    glEnableVertexAttribArray(0);
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, 0);
    glBindBuffer(GL_ARRAY_BUFFER, pbo_vel);
    glEnableVertexAttribArray(1);
    glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, 0);
    glBindBuffer(GL_ARRAY_BUFFER, 0);

    // Index buffer + CPU scratch for the tile-coherent draw order (rebuild_sort_order).
    glGenBuffers(1, &ibo);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ibo);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, (size_t)PBwidth * PBheight * sizeof(unsigned int), NULL, GL_DYNAMIC_DRAW);
    {
        int tilesW = (width + sort_tile - 1) / sort_tile;
        int tilesH = (height + sort_tile - 1) / sort_tile;
        sort_idx = (unsigned int *)malloc((size_t)P * sizeof(unsigned int));
        sort_key = (int *)malloc((size_t)P * sizeof(int));
        sort_counts = (int *)malloc(((size_t)tilesW * tilesH + 1) * sizeof(int));
        // Identity order until the first rebuild, so frame 0 still draws every particle.
        for (int i = 0; i < P; i++)
            sort_idx[i] = (unsigned int)i;
        glBufferSubData(GL_ELEMENT_ARRAY_BUFFER, 0, (size_t)P * sizeof(unsigned int), sort_idx);
    }

    // Core's own slots. The pool starts empty on every buffer_setup, including a --core.gl-refresh
    // context recreate, where the old names died with the old context.
    pool_reset();
    positionBufferIndex1 = pool_claim(1);
    positionBufferIndex2 = pool_claim(1);
    velocityBufferIndex1 = pool_claim(1);
    velocityBufferIndex2 = pool_claim(1);
    renderBufferIndex = pool_claim(1);
    repelBufferIndex = pool_claim(1);

    // Repel buffer (r = force potential, g = dispersion). Single-buffered and
    // stateless: fully cleared + re-splatted every sim_step (see the repel pass), so it
    // has no "other" half and never flips. Its own resolution knob (--core.repel-downsample), decoupled
    // from the density buffer -- this sets the hardness floor for the edge/exclusion/mouse repel.
    repelBuffer.minmag_filter = GL_LINEAR;
    repelBuffer.wrap_st = GL_REPEAT; // matches density/trail: velocity.frag samples it the same way
    repelBuffer.dim = BE_4D;         // r potential, g dispersion, ba drag vector
    buffer_allocate(&repelBuffer, current, repelBufferIndex, repel_width, repel_height, NULL);

    positionBuffer.minmag_filter = GL_NEAREST;
    positionBuffer.dim = BE_2D;
    velocityBuffer.minmag_filter = GL_NEAREST;
    velocityBuffer.dim = BE_2D;

    // Restore path (--core.gl-refresh): if a snapshot is queued and its dims match, upload
    // the saved particle/trail state instead of generating a fresh field, so the sim
    // continues across a GL-context recreate. NOTE: mismatched dims (only possible if
    // the surface resized between snapshot and setup) fall through to a fresh field.
    // Only the shared substrate's dims are core's business; each sim validates its own blob.
    bool restoring = g_restore && g_restore->pb_w == PBwidth && g_restore->pb_h == PBheight;

    // Position buffer 1.
    // --core.seed makes the initial positions (and so the whole deterministic gpu evolution)
    // reproducible, which is what lets before/after dumps actually be comparable.
    const float *pos_data;
    float *positions = NULL;
    if (restoring) {
        pos_data = g_restore->pos;
    } else {
        srand(rng_seed ? (unsigned int)rng_seed : (unsigned int)time(NULL));
        float margin = (sim_width < sim_height ? sim_width : sim_height) * 0.1f;
        float noise_seed = rng_seed ? (float)rng_seed : (10 * get_time() + frand01());
        fprintf(stdout, "Generating starting positions...\n");
        int N = PBwidth * PBheight * 2;
        positions = (float *)malloc(N * sizeof(float));
        for (int i = 0; i < N; i += 2) {
            // Generate random position in unit range
            vec2 position = {frand01(), frand01()};
            vec3 p0 = {10 * position[0], 10 * position[1], noise_seed};
            vec3 p1 = {10 * position[0], 10 * position[1], noise_seed + 1};
            vec2 noise = {glm_perlin_vec3(p0), glm_perlin_vec3(p1)};
            position[0] += (float)init_warp * fmodf(noise[0], 1.0f);
            position[1] += (float)init_warp * fmodf(noise[1], 1.0f);

            // Cull position to the center of the screen
            position[0] = position[0] * (sim_width - 2 * margin) + margin;
            position[1] = position[1] * (sim_height - 2 * margin) + margin;

            positions[i] = position[0];
            positions[i + 1] = position[1];
        }
        pos_data = positions;
    }

    fprintf(stdout, "Allocating buffers...\n");
    buffer_allocate(&positionBuffer, current, positionBufferIndex1, PBwidth, PBheight, (const char *)pos_data);

    free(positions);

    // Position buffer 2, velocity buffers 1 and 2. On restore, velocity
    // current gets the saved velocities (else zero); the "other" halves stay zero (the
    // double buffer flips fill them within a frame).
    buffer_allocate(&positionBuffer, other, positionBufferIndex2, PBwidth, PBheight, NULL);
    buffer_allocate(&velocityBuffer, current, velocityBufferIndex1, PBwidth, PBheight,
                    restoring ? (const char *)g_restore->vel : NULL);
    buffer_allocate(&velocityBuffer, other, velocityBufferIndex2, PBwidth, PBheight, NULL);

    // Internal render target (rgba colour) at render resolution.
    // GL_NEAREST so the upscale to the window stays crisp/pixelated at low res.
    renderBuffer.minmag_filter = GL_NEAREST;
    renderBuffer.wrap_st = GL_CLAMP_TO_EDGE;
    renderBuffer.dim = BE_4D;
    buffer_allocate(&renderBuffer, current, renderBufferIndex, width, height, NULL);

    // The screen buffer is the window itself (the upscale target), in device px.
    screenBuffer.width = window_width;
    screenBuffer.height = window_height;

    // The active sim's own fields, after the shared texture/framebuffer pool exists. It reads
    // the snapshot for its own restore, so this must run before the pointer is consumed below.
    for (size_t i = 0; i < NUM_SIMS; i++) {
        const SimBlob *blob = (restoring && i < MAX_SIMS) ? &g_restore->sims[i] : NULL;
        sims[i]->buffers(blob ? blob->data : NULL, blob ? blob->size : 0);
    }
    if (restoring && g_restore->active >= 0 && (size_t)g_restore->active < NUM_SIMS)
        active_sim = sims[g_restore->active]; // a refresh must not silently switch sims

    g_restore = NULL; // one restore per snapshot
}

SimSnapshot *sim_snapshot(void) {
    size_t pv = (size_t)PB_width * PB_height * 2;
    SimSnapshot *s = (SimSnapshot *)calloc(1, sizeof *s);
    if (s) {
        s->pos = malloc(pv * sizeof(float));
        s->vel = malloc(pv * sizeof(float));
    }
    if (!s || !s->pos || !s->vel) {
        fprintf(stderr, "goo: out of memory for sim snapshot\n");
        exit(EXIT_FAILURE);
    }
    s->pb_w = PB_width;
    s->pb_h = PB_height;

    // The shared substrate: read the *current* halves of each double buffer (what sim_step
    // just wrote). Fields derived from these -- density, and anything rebuilt every few frames
    // -- are deliberately not saved; they regenerate within a frame or two of the restore.
    glBindFramebuffer(GL_READ_FRAMEBUFFER, pool_framebuffers.items[positionBuffer.current]);
    glReadPixels(0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, s->pos);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, pool_framebuffers.items[velocityBuffer.current]);
    glReadPixels(0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, s->vel);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);

    // Every resident sim's own state, not just the active one's -- a switch after the restore
    // must find the other sim as it was left, not re-rolled.
    for (size_t i = 0; i < NUM_SIMS && i < MAX_SIMS; i++) {
        if (sims[i] == active_sim)
            s->active = (int)i;
        if (sims[i]->snapshot)
            s->sims[i].data = sims[i]->snapshot(&s->sims[i].size);
    }
    return s;
}

void sim_restore(const SimSnapshot *s) { g_restore = s; }

void sim_snapshot_free(SimSnapshot *s) {
    if (!s)
        return;
    free(s->pos);
    free(s->vel);
    for (size_t i = 0; i < MAX_SIMS; i++)
        free(s->sims[i].data);
    free(s);
}

void sim_teardown_cpu(void) {
    free(sort_idx);
    sort_idx = NULL;
    free(sort_key);
    sort_key = NULL;
    free(sort_counts);
    sort_counts = NULL;
    for (size_t i = 0; i < NUM_SIMS; i++) {
        // Every sim's exposure scratch, not just the active one's -- they are all resident.
        free(sims[i]->density->headroom.readback);
        sims[i]->density->headroom.readback = NULL;
        if (sims[i]->teardown_cpu)
            sims[i]->teardown_cpu(); // optional: a sim with no CPU scratch of its own leaves it NULL
    }
}

//============================================================
//
//   ####  ##  ###    ###      ###    #####  ##
//  ##     ##  ## #  # ##      ## ##   ##  ##  ##
//   ###   ##  ##  ##  ##     ##   ##  #####   ##
//     ##  ##  ##      ##     #######  ##      ##
//  ####   ##  ##      ##     ##   ##  ##      ##
//
//============================================================

void sim_set_dims(int logical_w, int logical_h, int fb_w, int fb_h) {
    // The actual framebuffer (window_width/height) is in device pixels and is only
    // the final upscale target. The sim renders at render resolution (width/height)
    // = logical points / render_scale, so on hidpi/retina it renders at 1x rather
    // than paying the device-pixel fill cost. The mouse comes in logical points, so
    // remember the ratio to scale it into render space.
    window_width = fb_w;
    window_height = fb_h;
    width = (int)((float)logical_w / render_scale + 0.5f);
    height = (int)((float)logical_h / render_scale + 0.5f);
    // Sim space is logical/SIM_SCALE, decoupled from the render resolution above.
    sim_width = (int)((double)logical_w / SIM_SCALE + 0.5);
    sim_height = (int)((double)logical_h / SIM_SCALE + 0.5);
    mouse_scale = (logical_w > 0) ? (float)sim_width / (float)logical_w : 1.0f;
}

// On resize the sim coordinate space (sim_width/height) changes, but the particle
// positions stored in positionBuffer don't -- left as-is they stay clumped in a corner
// of the new (usually larger) field and only slowly diffuse out to fill it. Scale every
// stored position by the size ratio so the field keeps filling the window immediately.
// Readback + reupload on the cpu: resizes are rare, so the cost is irrelevant, and it
// avoids standing up a dedicated scale shader pass. velocities are left unscaled -- they
// re-equilibrate within a few frames and the visible artefact is purely positional.
static void rescale_positions(float kx, float ky) {
    int n = PB_width * PB_height;
    float *p = (float *)malloc((size_t)n * 2 * sizeof(float));
    if (!p)
        return;
    glBindFramebuffer(GL_READ_FRAMEBUFFER, positionBuffer.framebuffers->items[positionBuffer.current]);
    glReadPixels(0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, p);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
    for (int i = 0; i < n; i++) {
        p[i * 2] *= kx;
        p[i * 2 + 1] *= ky;
    }
    // keep textures[i] bound to unit i (the convention buffer_allocate sets up).
    glActiveTexture(GL_TEXTURE0 + positionBuffer.current);
    glBindTexture(GL_TEXTURE_2D, positionBuffer.textures->items[positionBuffer.current]);
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, p);
    free(p);
}

void sim_resize(int logical_w, int logical_h, int fb_w, int fb_h) {
    int old_sim_w = sim_width, old_sim_h = sim_height;
    sim_set_dims(logical_w, logical_h, fb_w, fb_h);
    if (old_sim_w > 0 && old_sim_h > 0 && (sim_width != old_sim_w || sim_height != old_sim_h))
        rescale_positions((float)sim_width / (float)old_sim_w, (float)sim_height / (float)old_sim_h);
    glViewport(0, 0, window_width, window_height);
    updateShaderWindowShape(width, height);
    buffer_reallocate(&renderBuffer, current, width, height);
    // Repel buffer has its own resolution knob (--core.repel-downsample); recompute + realloc.
    repel_width = width / repelBufferDownsampling + 1;
    repel_height = height / repelBufferDownsampling + 1;
    buffer_reallocate(&repelBuffer, current, repel_width, repel_height);
    for (size_t i = 0; i < NUM_SIMS; i++)
        sims[i]->resize(); // each sim's own fields
    // sort_counts is sized by tile count, which grows with width/height -- resize it too,
    // else rebuild_sort_order writes past the end (sort_idx/sort_key are P-sized, stable).
    {
        int tilesW = (width + sort_tile - 1) / sort_tile;
        int tilesH = (height + sort_tile - 1) / sort_tile;
        sort_counts = (int *)realloc(sort_counts, ((size_t)tilesW * tilesH + 1) * sizeof(int));
    }
}

void sim_setup(void) {
    // Density is the default oracle for P: unless -p was given explicitly, derive P from the
    // logical (screen) area now that the front-end has set width/height. logical area =
    // width*height*render_scale^2 (width = logical/render_scale), so the count is independent of
    // both sim size and -r. Must run before buffer_setup, which consumes P.
    if (!p_given) {
        double logical_area = (double)width * height * render_scale * render_scale;
        P = (int)(init_density * logical_area + 0.5);
        if (P < 1)
            P = 1;
    }

    // Context state every pass below depends on: the point passes size their points from
    // the vertex shader, and the screen/debug passes composite SRC_ALPHA over. Set here
    // rather than left to the front-ends -- it's the core's own precondition, and a
    // front-end that omitted it rendered opaque with no error.
    glEnable(GL_PROGRAM_POINT_SIZE);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

    sim_build_pass_table();

    shader_setup();
    buffer_setup();
    // Static (per-run) uniforms, after shaders + buffers exist. The per-frame ones
    // (mouse, buffer indices, epoch) are set inside sim_step.
    for (size_t i = 0; i < NUM_SIMS; i++)
        sims[i]->static_uniforms();
}

// One physics + render frame. Runs the velocity / position / density / trail / screen
// passes into renderBuffer (render resolution). Does NOT present -- the front-end calls
// sim_present (or does its own readback) afterwards.
void sim_step(int epoch_counter, const float mouse_position[2], const float mouse_velocity[2]) {
    if (profile) {
        glFinish();
        pass_t0 = get_time();
    }

    // Interaction pass: fully clear + re-splat the force-potential (r) / dispersion (g) field
    // from scratch every frame -- bounding-rect edges, exclusion rects (--core.exclusions), and the
    // mouse. Single-buffered and stateless (no flip), so it's just a plain overwrite. Must run
    // BEFORE the velocity pass, which reads the fully-written result for this frame.
    // No blending: this pass fully clears and overwrites the field every frame, so the
    // fragment's alpha is DATA (the drag vector's y), not coverage. Leaving SRC_ALPHA on
    // multiplies the whole field by it and the walls vanish whenever the cursor is still.
    glDisable(GL_BLEND);
    shader_use(&repelShader);
    shader_set_uniform_vec(&repelShader, "mouse_position", 2, (float *)mouse_position);
    shader_set_uniform_vec(&repelShader, "mouse_velocity", 2, (float *)mouse_velocity);
    shader_set_uniform_vec(&repelShader, "buffer_shape", 2, (float[]){(float)repel_width, (float)repel_height});
    shader_set_uniform_float(&repelShader, "mouse_drag", (float)mouseDrag);
    buffer_bind(&repelBuffer, current);
    buffer_update(&repelBuffer);
    glEnable(GL_BLEND);

    // Which sim's physics runs this frame. Outside a transition it is simply the active one;
    // inside, error diffusion picks whole frames so the incoming sim's share matches the ramp.
    const SimDef *stepping = active_sim;
    float t = 1.0f; // incoming sim's share, 0 -> 1 across the transition
    if (transition_from) {
        t = switchFrames > 0 ? (float)transition_frame / (float)switchFrames : 1.0f;
        transition_acc += t;
        if (transition_acc >= 1.0)
            transition_acc -= 1.0; // this frame belongs to the incoming sim
        else
            stepping = transition_from;
    }

    // The sim advances its own per-particle state (velocity, position).
    stepping->integrate(epoch_counter, mouse_position, mouse_velocity);

    // Copy the fresh position + velocity textures into PBOs (FBO -> GL_PIXEL_PACK_BUFFER,
    // stays on-gpu). Velocity into pbo_vel; position into this frame's ring slot (async).
    glBindFramebuffer(GL_READ_FRAMEBUFFER, positionBuffer.framebuffers->items[positionBuffer.current]);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, pbo_pos[pos_pbo_ring]);
    glReadPixels(0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, 0);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, velocityBuffer.framebuffers->items[velocityBuffer.current]);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, pbo_vel);
    glReadPixels(0, 0, PB_width, PB_height, GL_RG, GL_FLOAT, 0);
    glBindBuffer(GL_PIXEL_PACK_BUFFER, 0);
    glBindFramebuffer(GL_READ_FRAMEBUFFER, 0);
    // Re-point the position vertex attribute (loc 0) at the slot just filled, so the point
    // passes below draw THIS frame's positions. (The VAO is still bound from buffer_setup.)
    glBindBuffer(GL_ARRAY_BUFFER, pbo_pos[pos_pbo_ring]);
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, 0);
    glBindBuffer(GL_ARRAY_BUFFER, 0);
    // Fence this slot's readback so rebuild_sort_order can poll (a few frames later) whether it
    // has landed before mapping it unsynchronized. Replace the slot's previous fence.
    if (pos_fence[pos_pbo_ring])
        glDeleteSync(pos_fence[pos_pbo_ring]);
    pos_fence[pos_pbo_ring] = glFenceSync(GL_SYNC_GPU_COMMANDS_COMPLETE, 0);
    LAP(0); // rdbk

    // Rebuild the tile-coherent draw order. rebuild_sort_order maps the OLDEST ring slot
    // unsynchronized, gated on that slot's fence -- it self-skips until the ring has a landed
    // readback (startup) or if the readback is momentarily behind, keeping the prior order.
    if (sort_every > 0 && epoch_counter % sort_every == 0) {
        double s0 = profile ? get_time() : 0; // CPU-side, so no glFinish -- but still --window.profile only
        rebuild_sort_order();
        if (profile)
            pass_ms[1] = ema(pass_ms[1], (get_time() - s0) * 1000.0, 0.1); // sort
    }
    pos_pbo_ring = (pos_pbo_ring + 1) % POS_PBO_N;

    // The sim splats/decays its own field buffers. Runs after the block above, which
    // repointed vertex attribute 0 at the positions integrate() just produced.
    stepping->fields(epoch_counter);

    // Screen rendering pass -> internal render target (render resolution)
    buffer_bind(&renderBuffer, current);
    if (transition_from) {
        // Both sims draw, weighted -- see the note by the transition state. The outgoing sim
        // first so the incoming one composites over it as its share grows.
        transition_from->draw(1.0f - t);
        active_sim->draw(t);
    } else {
        active_sim->draw(1.0f);
    }

    // Debug overlay: composite flat translucent-green fills into the render target (so they show
    // up in --window.dump/--window.headless which read renderBuffer, and regardless of particle density). Gated
    // on the C-side flags -> a normal run pays nothing. renderBuffer is still bound and blend is on
    // (SRC_ALPHA over). Fullscreen quad via gl_VertexID; the bound VAO needs no attributes.
    //  - exclusions_debug: fill each exclusion rect (rects live on the shader from push_exclusion_rects).
    //  - mouse_debug: fill the mouse repel-interaction disc. Radius is the same speed-scaled clamp the
    //    physics uses (repel.frag: clamp(speed*3, 30, 75) sim px -- constants mirrored here for the
    //    viz). sim px -> render px is * SIM_SCALE/render_scale; sim y is top-origin, gl_FragCoord is
    //    bottom-origin, so y flips. A parked mouse sits far off-screen, so its disc covers nothing.
    //  - edge_debug: fill the bounding-edge repel bands (render_shape lets the shader size them).
    //  - mouse_debug also draws the opaque cursor dot + a motion line from last frame's cursor
    //    position (mouse_prev = current - per-frame velocity) to the current one.
    if (exclusions_debug || mouse_debug || edge_debug) {
        shader_use(&debugShader);
        shader_set_uniform_int(&debugShader, "draw_exclusions", exclusions_debug ? 1 : 0);
        shader_set_uniform_int(&debugShader, "draw_mouse", mouse_debug ? 1 : 0);
        shader_set_uniform_int(&debugShader, "draw_edges", edge_debug ? 1 : 0);
        shader_set_uniform_vec(&debugShader, "render_shape", 2, (float[]){(float)width, (float)height});
        // repel buffer's own resolution -- lets the overlay snap to that (downsampled) texel
        // grid, so it shows the ACTUAL coarse field the physics samples, not the ideal shape.
        shader_set_uniform_vec(&debugShader, "buffer_shape", 2, (float[]){(float)repel_width, (float)repel_height});
        shader_set_uniform_int(&debugShader, "repel_buffer", repelBuffer.current);
        float mspeed = sqrtf(mouse_velocity[0] * mouse_velocity[0] + mouse_velocity[1] * mouse_velocity[1]);
        // ease-out radius, mirrors repel.frag/velocity.frag: fast rise at low speed, plateau at max.
        float mt = fminf(mspeed / 20.0f, 1.0f);
        float mrad_sim = 30.0f + (75.0f - 30.0f) * (1.0f - (1.0f - mt) * (1.0f - mt));
        float k = (float)(SIM_SCALE / render_scale);
        shader_set_uniform_vec(&debugShader, "mouse_center", 2,
                               (float[]){mouse_position[0] * k, (float)height - mouse_position[1] * k});
        shader_set_uniform_vec(&debugShader, "mouse_prev", 2,
                               (float[]){(mouse_position[0] - mouse_velocity[0]) * k,
                                         (float)height - (mouse_position[1] - mouse_velocity[1]) * k});
        shader_set_uniform_float(&debugShader, "mouse_radius", mrad_sim * k);
        glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
    }
    LAP(2); // draw (the sim's screen pass + the debug overlay)

    if (transition_from) {
        // Drop anything requested mid-transition rather than letting it sit: an un-cleared flag
        // would fire on the very next frame, starting the reverse transition with zero settled
        // frames on the sim that just arrived.
        switch_requested = 0;
        if (++transition_frame > switchFrames) {
            fprintf(stdout, "sim: now %s\n", active_sim->name);
            transition_from = NULL;
        }
    } else {
        // --core.switch-every: the timer shares the trigger path with tab and SIGUSR1, and like
        // them cannot fire mid-transition. Measured from the end of the last one, not from start,
        // so the cadence is time SETTLED rather than time between starts.
        static double last_switch = 0;
        if (last_switch == 0)
            last_switch = get_time();
        if (switchEvery > 0 && get_time() - last_switch >= switchEvery)
            switch_requested = 1;
        if (switch_requested) {
            switch_requested = 0;
            last_switch = get_time();
            sim_begin_switch();
        }
    }
}

void sim_present(void) {
    // Timed on its own mark rather than through LAP: present runs outside sim_step, so
    // the frame's lap chain has already ended by the time the front-end calls this.
    double t0 = 0;
    if (profile) {
        glFinish();
        t0 = get_time();
    }

    // Upscale pass: nearest-sample the render target onto the window. Blending
    // is off -- the render target already holds particles composited over black,
    // so this is a straight opaque blit (chunky when render_scale > 1).
    buffer_bind(&screenBuffer, screen);
    glDisable(GL_BLEND);
    shader_use(&upscaleShader);
    shader_set_uniform_int(&upscaleShader, "source_buffer", renderBuffer.current);
    // Source crop: identity (full field) for every normal present; macwp -m=-2 sets a
    // centred sub-rect so each monitor shows its own 1:1 slice of the shared field.
    shader_set_uniform_vec(&upscaleShader, "uv_min", 2, present_uv_min);
    shader_set_uniform_vec(&upscaleShader, "uv_max", 2, present_uv_max);
    buffer_update(&screenBuffer);
    glEnable(GL_BLEND);

    if (profile) {
        glFinish();
        pass_ms[3] = ema(pass_ms[3], (get_time() - t0) * 1000.0, 0.1); // up
    }
}