local _tline3d = tline3d

-- simple deterministic RNG for shake/glitch
local _rng = 0x12345678
local function irnd(n)
  _rng = (_rng * 1103515245 + 12345) & 0x7fffffff
  return _rng % n
end

local CFG = {
  w = 240, h = 135, vid = 3,

  bg_col   = 0,
  hud_col  = 63,

  focal = 192,
  near  = 0.15,
  cam_roll_speed = 0.45,  -- radians per second
  bg_rot_speed   = -0.12,  -- background rotation speed (radians per second)

  -- wall geometry
  grid_x    = 4.4,
  grid_y0   = -4.8,
  grid_y1   =  4.8,
  grid_z0   =  7.0,
  grid_z1   = 18.7,

  -- wall iso-surface
  wall_str   = 0.5,
  wall_thick = 0.2,

  -- blob motion
  num_blobs  = 4,
  blob_cycles = {4.2, 5.8, 7.5, 9.0,12.0, 23.5},  -- per-blob traverse time
  attach_t   = 0.18, -- time spent fading in/out at start/end of path
  phase_skew = 0.001, -- offset between path progress and blob position (0=no skew, 1=full path delay)

  r_min = 1.00, -- minimum blob radius
  r_max = 2.50, -- maximum blob radius

  bend_y = 1.35, -- max vertical control point offset (blob path curvature)
  bend_z = 1.55, -- max depth control point offset (blob path curvature)
  wob_y  = 0.88, -- max vertical wobble amplitude
  wob_z  = 0.84, -- max depth wobble amplitude 

  -- scalar field
  iso  = 0.01, -- minimum field value to be considered "inside" the surface (iso-surface threshold)
  soft = 0.2, -- softening factor added to squared distance when evaluating field strength (prevents singularities and allows smoother blending)

  -- wall anchor cluster
  anchor_gain     = 1.10, -- multiplier on anchor field strength (relative to blob center) to make them more/less influential
  anchor_side_w   = 0.46, -- relative strength of side anchors (0=only center, 1=same as center)
  anchor_spread_y = 0.60, -- vertical spread of side anchors as multiple of blob radius
  anchor_spread_z = 0.78, -- depth spread of side anchors as multiple of blob radius
  anchor_inset    = 0.42, -- how far from wall center the anchor cluster is placed (relative to grid_x)

  -- voxel lattice (finer grid)
  nx = 24,
  ny = 24,
  nz = 14,

  x0 = -5.8, x1 =  5.8,
  y0 = -5.8, y1 =  5.8,
  z0 =  5.5, z1 = 20.8,

  vox_scale = 0.35,
}

local sin  = math.sin
local cos  = math.cos
local abs  = math.abs
local min  = math.min
local max  = math.max
local flr  = math.floor
local sqrt = math.sqrt

local cx, cy
local NX, NY, NZ
local DX, DY, DZ
local VOX

local TOTAL

local blob_defs = {}   -- path definition per blob
local blob_ids  = {}   -- current path id per blob (avoid recalc)
local blobs     = {}   -- evaluated blob state per frame

local active_iy0, active_iy1 = 0, 0
local active_iz0, active_iz1 = 0, 0
local active_ix0, active_ix1 = 0, 0

local sxs, sys, szs, sss = {}, {}, {}, {}
local src_n = 0

-- pre-computed coordinate arrays (Lua tables for draw loop)
local vx_arr = {}
local vy_arr = {}
local vz_arr = {}

-- userdata arrays for batch field evaluation
local cell_x_ud, cell_y_ud, cell_z_ud
local wall_field_ud
local blob_mask_ud  -- 1.0 between walls, 0.0 outside
local field_ud
local tmp1_ud, tmp2_ud
local r2_base_ud  -- pre-computed x²+y²+z²+soft (constant)

-- depth LUT: maps iz (0..NZ-1) to a depth brightness 0..1
local depth_lut = {}

-- base palette colors (stored by setup_palette for fade)
local base_pal = {}

-- HSV helpers for hue rotation
local function rgb_to_hsv(r, g, b)
  r, g, b = r / 255, g / 255, b / 255
  local mx = max(r, max(g, b))
  local mn = min(r, min(g, b))
  local d = mx - mn
  local h, s, v = 0, 0, mx
  if mx > 0 then s = d / mx end
  if d > 0 then
    if mx == r then     h = (g - b) / d % 6
    elseif mx == g then h = (b - r) / d + 2
    else                h = (r - g) / d + 4
    end
    h = h / 6
  end
  return h, s, v
end

local function hsv_to_rgb(h, s, v)
  h = h % 1.0
  if h < 0 then h = h + 1 end
  local c = v * s
  local x = c * (1 - abs(h * 6 % 2 - 1))
  local m = v - c
  local r, g, b = m, m, m
  local hi = flr(h * 6)
  if     hi == 0 then r = r + c; g = g + x
  elseif hi == 1 then r = r + x; g = g + c
  elseif hi == 2 then g = g + c; b = b + x
  elseif hi == 3 then g = g + x; b = b + c
  elseif hi == 4 then r = r + x; b = b + c
  else                r = r + c; b = b + x
  end
  return flr(min(255, r * 255)), flr(min(255, g * 255)), flr(min(255, b * 255))
end

-----------------------------
-- PALETTE
-----------------------------
local pal_entries = {}  -- loaded from pod.pal file
local bg_tex, bg_tw, bg_th, bg_iw, bg_ih  -- background texture (pow2) for tline3d

local function load_pal_file()
  local raw = fetch("metaplanes.pod.pal")
  if not raw then return end
  local txt = raw
  if type(raw) == "table" then txt = table.concat(raw, "\n") end
  if type(txt) ~= "string" then return end
  pal_entries = {}
  for token in txt:gmatch("[#]?([0-9a-fA-F][0-9a-fA-F][0-9a-fA-F][0-9a-fA-F][0-9a-fA-F][0-9a-fA-F])") do
    local v = tonumber(token, 16)
    if v then pal_entries[#pal_entries + 1] = v end
  end
end

local function setup_palette()
  for i = 0, 63 do
    local c = pal_entries[i + 1] or 0x050510
    pal(i, c, 2)
    base_pal[i] = c
  end
end

local function clamp(x,a,b)
  if x < a then return a end
  if x > b then return b end
  return x
end

local function lerp(a,b,t)
  return a + (b-a) * t
end

local function smooth(t)
  t = clamp(t, 0, 1)
  return t*t*(3-2*t)
end

local function frac(x)
  return x - flr(x)
end

local function hash(n)
  return frac(sin(n * 12.9898 + 78.233) * 43758.5453123)
end

local function randf(seed, a, b)
  return a + (b-a) * hash(seed)
end

local function bez(a,b,c,d,t)
  local it = 1 - t
  return it*it*it*a
       + 3*it*it*t*b
       + 3*it*t*t*c
       + t*t*t*d
end



local function proj(x,y,z)
  if z <= CFG.near then return nil end
  local s = CFG.focal / z
  return cx + x * s, cy + y * s, s
end

local function field_strength_from_r(r)
  return CFG.iso * (r*r + CFG.soft)
end

local function reset_sources()
  src_n = 0
end

local function add_source(x,y,z,s)
  src_n += 1
  sxs[src_n] = x
  sys[src_n] = y
  szs[src_n] = z
  sss[src_n] = s
end

local function add_anchor_cluster(side, y, z, r, w)
  local xc = side * CFG.grid_x
  local xi = xc - side * CFG.anchor_inset
  local oy = r * CFG.anchor_spread_y
  local oz = r * CFG.anchor_spread_z
  local s  = field_strength_from_r(r) * w * CFG.anchor_gain

  add_source(xc, y,      z,      s * 1.15)
  add_source(xi, y,      z,      s * 0.95)
  add_source(xc, y + oy, z,      s * CFG.anchor_side_w)
  add_source(xc, y - oy, z,      s * CFG.anchor_side_w)
  add_source(xc, y,      z + oz, s * CFG.anchor_side_w)
  add_source(xc, y,      z - oz, s * CFG.anchor_side_w)
end



local function build_blob_def(bi, id)
  if blob_ids[bi] == id then return end

  local prev = blob_defs[bi]  -- previous path (nil on first call)
  blob_ids[bi] = id

  local p = prev or {}  -- reuse existing table

  -- restrict blobs to the far third of the z range, near center
  local blob_z0 = CFG.grid_z0 + (CFG.grid_z1 - CFG.grid_z0) * (2/3)
  local blob_z1 = CFG.grid_z1
  local blob_zc = (blob_z0 + blob_z1) * 0.5  -- center of far third
  local z_off = 1.2  -- small offset from center

  if prev then
    -- chain: start where the previous path ended
    p.side = -prev.side        -- arrived at opposite wall, now depart from there
    p.sy   = prev.ey
    p.sz   = prev.ez
  else
    -- first path: random start near center of planes
    p.side = (hash(id * 71 + bi * 37) > 0.5) and 1 or -1
    p.sy = randf(id*101 + bi*17, CFG.grid_y0*0.68, CFG.grid_y1*0.68)
    p.sz = randf(id*131 + bi*23, blob_zc - z_off, blob_zc + z_off)
  end

  p.ey = randf(id*151 + bi*29, CFG.grid_y0*0.68, CFG.grid_y1*0.68)
  p.ez = randf(id*181 + bi*31, blob_zc - z_off, blob_zc + z_off)

  p.c1y = p.sy + randf(id*199 + bi*37, -CFG.bend_y, CFG.bend_y)
  p.c2y = p.ey + randf(id*223 + bi*41, -CFG.bend_y, CFG.bend_y)

  p.c1z = lerp(p.sz, p.ez, 0.33) + randf(id*239 + bi*43, -CFG.bend_z, CFG.bend_z)
  p.c2z = lerp(p.sz, p.ez, 0.66) + randf(id*251 + bi*47, -CFG.bend_z, CFG.bend_z)

  p.r = randf(id*271 + bi*53, CFG.r_min, CFG.r_max)
  p.skew = randf(id*281 + bi*59, -CFG.phase_skew, CFG.phase_skew)

  p.wob_ay = randf(id*307 + bi*61, 0.0, CFG.wob_y)
  p.wob_az = randf(id*331 + bi*67, 0.0, CFG.wob_z)
  p.wob_f  = randf(id*353 + bi*71, 0.85, 1.45)
  p.wob_p  = randf(id*379 + bi*73, 0.0, 6.28318)

  blob_defs[bi] = p
end

local function eval_blob(def, ph, t, out)
  local u = clamp(ph + def.skew, 0, 1)
  local e = smooth(u)

  local gx = CFG.grid_x
  local raw_x = lerp(def.side * gx, -def.side * gx, e)

  -- clamp blob center to stay inside walls (inset by wall_thick)
  local wall_in = CFG.wall_thick
  local x = clamp(raw_x, -gx + wall_in, gx - wall_in)

  local y = bez(def.sy, def.c1y, def.c2y, def.ey, e)
          + sin(t * def.wob_f + def.wob_p) * def.wob_ay * sin(e * 3.14159)

  local z = bez(def.sz, def.c1z, def.c2z, def.ez, e)
          + cos(t * (def.wob_f * 1.17) + def.wob_p * 1.53) * def.wob_az * sin(e * 3.14159)

  local src_att = 1 - smooth(clamp(u / CFG.attach_t, 0, 1))
  local dst_att = 1 - smooth(clamp((1-u) / CFG.attach_t, 0, 1))

  local r = def.r * (1 + 0.07 * sin(t * 1.67 + def.wob_p))

  -- fade based on how much of the blob's outer edge is past the wall
  -- blob is fully visible when its inner edge (center-r toward center) is far from wall
  -- blob is fully absorbed when its inner edge reaches the wall plane
  local dx_wall = gx - abs(x)         -- center distance to wall
  local edge_dist = dx_wall - r       -- inner edge distance (negative = edge past wall)
  local fade = clamp(edge_dist / r + 1, 0, 1)  -- 0 when edge at wall, 1 when one radius away

  out.side = def.side
  out.x = x; out.y = y; out.z = z
  out.r = r; out.fade = fade
  out.sy = def.sy; out.sz = def.sz
  out.ey = def.ey; out.ez = def.ez
  out.src_att = src_att
  out.dst_att = dst_att
end

local function build_field(t)
  reset_sources()

  local nb = CFG.num_blobs
  local cycles = CFG.blob_cycles

  for bi=1,nb do
    local cyc_i = t / cycles[bi]
    local id_i = flr(cyc_i)
    local ph_i = cyc_i - id_i

    build_blob_def(bi, id_i)
    if not blobs[bi] then blobs[bi] = {} end
    eval_blob(blob_defs[bi], ph_i, t, blobs[bi])

    local b = blobs[bi]
    add_source(b.x, b.y, b.z, field_strength_from_r(b.r) * b.fade)

    if b.src_att > 0.001 and b.fade < 0.95 then
      add_anchor_cluster(b.side,  b.sy, b.sz, b.r, b.src_att * (1 - b.fade))
    end
    if b.dst_att > 0.001 and b.fade < 0.95 then
      add_anchor_cluster(-b.side, b.ey, b.ez, b.r, b.dst_att * (1 - b.fade))
    end
  end

  -- AABB: compute active grid index ranges from source positions
  local eff_r = 3.0
  local bnd_x0, bnd_x1 = CFG.grid_x, -CFG.grid_x
  local bnd_y0, bnd_y1 = CFG.grid_y0, CFG.grid_y1
  local bnd_z0, bnd_z1 = CFG.grid_z0, CFG.grid_z1
  for si=1,src_n do
    local sxi, syi, szi = sxs[si], sys[si], szs[si]
    if sxi - eff_r < bnd_x0 then bnd_x0 = sxi - eff_r end
    if sxi + eff_r > bnd_x1 then bnd_x1 = sxi + eff_r end
    if syi - eff_r < bnd_y0 then bnd_y0 = syi - eff_r end
    if syi + eff_r > bnd_y1 then bnd_y1 = syi + eff_r end
    if szi - eff_r < bnd_z0 then bnd_z0 = szi - eff_r end
    if szi + eff_r > bnd_z1 then bnd_z1 = szi + eff_r end
  end
  -- walls always span full x range
  if bnd_x0 > -CFG.grid_x then bnd_x0 = -CFG.grid_x end
  if bnd_x1 <  CFG.grid_x then bnd_x1 =  CFG.grid_x end
  active_ix0 = max(0, flr((bnd_x0 - CFG.x0) / DX))
  active_ix1 = min(NX-1, flr((bnd_x1 - CFG.x0) / DX) + 1)
  active_iy0 = max(0, flr((bnd_y0 - CFG.y0) / DY))
  active_iy1 = min(NY-1, flr((bnd_y1 - CFG.y0) / DY) + 1)
  active_iz0 = max(0, flr((bnd_z0 - CFG.z0) / DZ))
  active_iz1 = min(NZ-1, flr((bnd_z1 - CFG.z0) / DZ) + 1)
end

-----------------------------
-- USERDATA-VECTORIZED FIELD EVALUATION + DRAW
-- All source contributions computed via batch userdata arithmetic
-- (no per-cell Lua source loop). Wall planes are part of the iso-surface.
-----------------------------

local function compute_and_draw_voxels_rocket(this)
  local _sxs, _sys, _szs, _sss = sxs, sys, szs, sss
  local _src_n = src_n
  local plane = NY * NX

  -- 1. Initialize blob field to zero
  field_ud:mul(0, true)

  -- 2. Batch evaluate sources with z-range culling
  -- Each source only processes z-planes within its effective radius,
  -- using offset+count parameters on in-place userdata ops.
  -- Contribution at z-distance dz: ss/(dz²+soft). Skip when < iso*0.1.
  local _t1 = tmp1_ud
  local _t2 = tmp2_ud
  local _r2b = r2_base_ud
  local z0_grid = CFG.z0
  local soft = CFG.soft
  local skip_th = CFG.iso * 0.1

  for si=1,_src_n do
    local sx, sy, sz, ss = _sxs[si], _sys[si], _szs[si], _sss[si]

    if ss >= 0.0001 then
      -- compute z-range where this source matters
      local max_dz2 = ss / skip_th - soft
      if max_dz2 > 0 then
        local max_dz = sqrt(max_dz2)
        local iz_lo = max(0, flr((sz - max_dz - z0_grid) / DZ))
        local iz_hi = min(NZ-1, flr((sz + max_dz - z0_grid) / DZ) + 1)
        if iz_hi >= iz_lo then
          local off = iz_lo * plane
          local cnt = (iz_hi - iz_lo + 1) * plane
          local src_r2 = sx*sx + sy*sy + sz*sz

          -- t1 = -2·sx·x
          _t1:copy(cell_x_ud, true, off, off, cnt)
          _t1:mul(-2*sx, true, 0, off, cnt)
          -- t2 = -2·sy·y
          _t2:copy(cell_y_ud, true, off, off, cnt)
          _t2:mul(-2*sy, true, 0, off, cnt)
          _t1:add(_t2, true, off, off, cnt)
          -- t2 = -2·sz·z
          _t2:copy(cell_z_ud, true, off, off, cnt)
          _t2:mul(-2*sz, true, 0, off, cnt)
          _t1:add(_t2, true, off, off, cnt)
          -- t1 += r2_base + src_r² → t1 = dist²+soft
          _t1:add(_r2b, true, off, off, cnt)
          _t1:add(src_r2, true, 0, off, cnt)

          -- accumulate: field += ss/t1 via (field·t1 + ss)/t1
          field_ud:mul(_t1, true, off, off, cnt)
          field_ud:add(ss, true, 0, off, cnt)
          field_ud:div(_t1, true, off, off, cnt)
        end
      end
    end
  end

  -- 3. Mask blob field, compute surf + blob_ratio via destination-operand trick
  -- Only process z-planes within the active draw range (offset+count)
  local post_off = active_iz0 * plane
  local post_cnt = (active_iz1 - active_iz0 + 1) * plane
  field_ud:mul(blob_mask_ud, true, post_off, post_off, post_cnt)  -- field = masked blob field
  local inv_iso = 1.0 / CFG.iso
  field_ud:add(wall_field_ud, tmp2_ud)            -- tmp2 = total (blob+wall), field unchanged
  tmp2_ud:add(0.001, tmp1_ud)                     -- tmp1 = total + eps
  field_ud:div(tmp1_ud, true)                     -- field = blob_ratio = blob/(total+eps)
  tmp2_ud:mul(inv_iso, tmp1_ud)                   -- tmp1 = surf = total * inv_iso

  -- 4. Draw back-to-front (painter's algorithm)
  local focal   = CFG.focal
  local vs      = VOX * CFG.vox_scale
  local near    = CFG.near
  local cam_ox  = this.cam_x
  local cam_oy  = this.cam_y
  local cam_oz  = this.cam_z
  local _cx, _cy = cx + cam_ox, cy + cam_oy
  local _vy, _vz = vy_arr, vz_arr
  local _surf   = tmp1_ud     -- pre-computed surf (for visibility + color)
  local _bratio = field_ud    -- pre-computed blob_ratio
  local _rf     = rectfill
  local _aiy0, _aiy1 = active_iy0, active_iy1
  local _aix0, _aix1 = active_ix0, active_ix1
  local _NX     = NX
  local x0_world = CFG.x0
  local SCR_W, SCR_H = CFG.w, CFG.h

  -- camera roll rotation
  local roll = this.t * CFG.cam_roll_speed
  local cr = cos(roll)
  local sr = sin(roll)

  -- screen-center brightness boost parameters
  local half_w, half_h = SCR_W * 0.5, SCR_H * 0.5
  local inv_max_r = 1.0 / sqrt(half_w * half_w + half_h * half_h)
  local _gx = CFG.grid_x
  local _inv_gx = 1.0 / _gx
  local _vx = vx_arr

  for iz=active_iz1, active_iz0, -1 do
    local z = _vz[iz] + cam_oz
    if z > near then
      local inv_z = focal / z
      local r = vs * inv_z
      local d62 = depth_lut[iz] * 62  -- pre-scaled depth*62
      local k_base = iz * plane

      -- adaptive stride: when projected rect overlaps neighbours, skip cells
      -- ratio = projected rect size / projected cell spacing
      local px_cell = DX * inv_z  -- projected cell spacing
      local ratio = r / px_cell
      local stride = 1
      if ratio > 1.5 then stride = 2 end
      if ratio > 3.0 then stride = 4 end
      -- scale rect to fill the gaps from skipped cells
      local r_draw = r * stride

      -- DDA: projected x advances by constant step per ix
      local px_step = px_cell * stride
      local px_start = (x0_world + _aix0 * DX) * inv_z
      -- screen coord increments per ix (roll-rotated)
      local sx_step = px_step * cr
      local sy_step = px_step * sr
      -- screen bounds for frustum cull (with margin for voxel radius)
      local scr_x0 = -r_draw
      local scr_x1 = SCR_W + r_draw
      local scr_y0 = -r_draw
      local scr_y1 = SCR_H + r_draw
      -- ix span count
      local ix_span = _aix1 - _aix0

      -- projected y cell spacing for iy stride
      local py_cell = DY * inv_z
      local y_ratio = r / py_cell
      local y_stride = 1
      if y_ratio > 1.5 then y_stride = 2 end
      if y_ratio > 3.0 then y_stride = 4 end

      for iy=_aiy0,_aiy1,y_stride do
        local row_off = k_base + iy * _NX
        -- hoist py projection + roll rotation out of ix loop
        local py = _vy[iy] * inv_z
        local sx_row = _cx + px_start * cr - py * sr
        local sy_row = _cy + px_start * sr + py * cr

        -- row-level frustum cull: check if entire row's sy range is off-screen
        local sy_end = sy_row + sy_step * ix_span
        local sy_lo = sy_row < sy_end and sy_row or sy_end
        local sy_hi = sy_row < sy_end and sy_end or sy_row
        if sy_hi >= scr_y0 and sy_lo <= scr_y1 then
          -- also check sx range for the row
          local sx_end = sx_row + sx_step * ix_span
          local sx_lo = sx_row < sx_end and sx_row or sx_end
          local sx_hi = sx_row < sx_end and sx_end or sx_row
          if sx_hi >= scr_x0 and sx_lo <= scr_x1 then

            for ix=_aix0,_aix1,stride do
              local sv = _surf[row_off + ix]
              if sv >= 1.0 then
                if sv > 2.0 then sv = 2.0 end
                local sa = sv - 1.0

                local br = _bratio[row_off + ix]
                if br > 1.0 then br = 1.0 end

                -- depth + center boost + blobs brighter away from walls
                local dx = sx_row - half_w
                local dy = sy_row - half_h
                local dist = sqrt(dx * dx + dy * dy) * inv_max_r
                local boost = 1.0 - dist * 0.5  -- center=1.0, corners=0.5
                -- wall proximity: 0 at wall, 1 at center
                local wall_frac = (_gx - abs(_vx[ix])) * _inv_gx
                if wall_frac < 0 then wall_frac = 0 elseif wall_frac > 1 then wall_frac = 1 end
                local col = flr(d62 * boost * (0.6 + 0.4 * br * wall_frac) + 1)
                if col < 1 then col = 1 elseif col > 63 then col = 63 end

                _rf(sx_row - r_draw, sy_row - r_draw, sx_row + r_draw, sy_row + r_draw, col)
              end
              sx_row = sx_row + sx_step
              sy_row = sy_row + sy_step
            end

          end -- sx range cull
        end -- sy range cull
      end
    end
  end
end

------------------------------------------------------------
-- rocket parameters (active values, set via setparams)
--
-- params[2] = time         (animation time, drives everything)*
-- params[3] = focal        (camera focal length, 192 = default)*
-- params[4] = cam_roll     (camera roll speed, 0.45 = default)*
-- params[5] = cam_x        (camera x offset, 0 = default)*
-- params[6] = cam_y        (camera y offset, 0 = default)*
-- params[7] = cam_z        (camera z offset, 0 = default)*
-- params[8] = vox_scale    (rectangle size multiplier, 0.35 = default)*
-- params[9] = fade         (0=black, 1=normal, 2=white)*
-- params[10]= shake        (screen shake amplitude, 0=off)*
-- params[11]= glitch       (glitch intensity, 0=off)*
-- params[12]= hue          (hue rotation 0..1, 0 = default)*
-- params[13]= contrast     (contrast 0..2, 1 = default, 0=flat grey, 2=harsh)*
------------------------------------------------------------

metaplanes = {
  t          = 0,
  focal      = 192,
  cam_roll   = 0.45,
  cam_x      = 0,
  cam_y      = 0,
  cam_z      = 0,
  vox_scale  = 0.35,
  fade       = 1.0,
  shake      = 0,
  glitch     = 0,
  hue        = 0,
  contrast   = 1.0,
  _pal_dirty = false,
}

metaplanes.load = function(this)
  vid(CFG.vid)
  load_pal_file()
  setup_palette()
  palt(0, false)

  -- load background image into pow2 texture for tline3d rotation
  local obj = fetch("metaplanes.pod")
  if obj then
    local raw = unpod(obj)
    if not raw and type(obj) == "table" then
      raw = obj.bmp or obj.image or obj.img
    end
    if raw then
      local iw, ih = raw:width(), raw:height()
      local pw, ph = 1, 1
      while pw < iw do pw *= 2 end
      while ph < ih do ph *= 2 end
      bg_tex = userdata("u8", pw, ph)
      set_draw_target(bg_tex)
      cls(0)
      palt(0, false)
      sspr(raw, 0, 0, iw, ih, 0, 0, iw, ih)
      set_draw_target()
      bg_tw, bg_th = pw, ph
      bg_iw, bg_ih = iw, ih
    end
  end

  cx = CFG.w * 0.5
  cy = CFG.h * 0.5

  NX, NY, NZ = CFG.nx, CFG.ny, CFG.nz
  TOTAL = NX * NY * NZ
  DX = (CFG.x1 - CFG.x0) / (NX - 1)
  DY = (CFG.y1 - CFG.y0) / (NY - 1)
  DZ = (CFG.z1 - CFG.z0) / (NZ - 1)
  VOX = min(DX, min(DY, DZ))

  -- Pre-compute coordinate arrays (Lua tables for draw loop)
  for ix=0,NX-1 do vx_arr[ix] = CFG.x0 + ix * DX end
  for iy=0,NY-1 do vy_arr[iy] = CFG.y0 + iy * DY end
  for iz=0,NZ-1 do vz_arr[iz] = CFG.z0 + iz * DZ end

  -- Pre-compute depth brightness LUT (close=bright, far=dark)
  for iz=0,NZ-1 do
    local z = vz_arr[iz]
    local t = (z - CFG.z0) / (CFG.z1 - CFG.z0)  -- 0=near, 1=far
    if t < 0 then t = 0 elseif t > 1 then t = 1 end
    depth_lut[iz] = 1.0 - t * 0.9  -- near=1.0, far=0.1
  end

  -- Allocate userdata arrays for batch field evaluation
  cell_x_ud     = userdata("f64", TOTAL)
  cell_y_ud     = userdata("f64", TOTAL)
  cell_z_ud     = userdata("f64", TOTAL)
  wall_field_ud = userdata("f64", TOTAL)
  blob_mask_ud  = userdata("f64", TOTAL)
  field_ud      = userdata("f64", TOTAL)
  tmp1_ud       = userdata("f64", TOTAL)
  tmp2_ud       = userdata("f64", TOTAL)
  r2_base_ud    = userdata("f64", TOTAL)

  -- Fill coordinate arrays and pre-compute wall field
  local gx     = CFG.grid_x
  local gy0    = CFG.grid_y0
  local gy1    = CFG.grid_y1
  local gz0    = CFG.grid_z0
  local gz1    = CFG.grid_z1
  local wstr   = CFG.wall_str
  local wthick = CFG.wall_thick

  local k = 0
  for iz=0,NZ-1 do
    local z = vz_arr[iz]
    for iy=0,NY-1 do
      local y = vy_arr[iy]
      for ix=0,NX-1 do
        local x = vx_arr[ix]
        cell_x_ud[k] = x
        cell_y_ud[k] = y
        cell_z_ud[k] = z
        r2_base_ud[k] = x*x + y*y + z*z + CFG.soft

        -- Wall field: thin plane at x = ±grid_x within wall y/z bounds
        local wf = 0
        if y >= gy0 and y <= gy1 and z >= gz0 and z <= gz1 then
          local dl = abs(x + gx)    -- distance to left wall (-grid_x)
          local dr = abs(x - gx)    -- distance to right wall (+grid_x)
          local d = min(dl, dr)
          if d < wthick then
            local t = 1 - d / wthick
            wf = wstr * t * t
          end
        end
        wall_field_ud[k] = wf

        -- blob mask: 1.0 between walls, fade to 0.0 past them
        local bm = 1.0
        local mask_margin = 0.6
        local ax = abs(x)
        if ax > gx then
          bm = clamp(1.0 - (ax - gx) / mask_margin, 0, 1)
        end
        blob_mask_ud[k] = bm

        k = k + 1
      end
    end
  end

  -- pre-warm: build all blob definitions and evaluate field at t=0
  -- this forces table allocations and path generation to happen at load time
  -- instead of causing a CPU spike on the first draw frame
  build_field(0)
end

metaplanes.init = function(this)
  setup_palette()
  this._pal_dirty = false
end

metaplanes.setparams = function(this, params)
  this.t          = params[2] or 0
  this.focal      = params[3] or 192
  this.cam_roll   = params[4] or 0.45
  this.cam_x      = params[5] or 0
  this.cam_y      = params[6] or 0
  this.cam_z      = params[7] or 0
  if params[8] == nil then this.vox_scale = 0.35 else this.vox_scale = params[8] end
  this.fade       = params[9] or 1.0
  this.shake      = params[10] or 0
  this.glitch     = params[11] or 0
  if params[12] == nil then this.hue = 0 else this.hue = params[12] end
  this.contrast   = params[13] or 1.0
end

metaplanes.update = function(this)
  local fade = this.fade
  local hue_shift = this.hue
  local con = this.contrast
  local needs_transform = fade ~= 1.0 or hue_shift ~= 0 or con ~= 1.0

  if needs_transform then
    for i = 0, 63 do
      local c = base_pal[i]
      if c then
        local r = (c >> 16) & 0xff
        local g = (c >> 8) & 0xff
        local b = c & 0xff

        -- hue rotation
        if hue_shift ~= 0 then
          local h, s, v = rgb_to_hsv(r, g, b)
          h = h + hue_shift
          r, g, b = hsv_to_rgb(h, s, v)
        end

        -- contrast (around mid-grey 128)
        if con ~= 1.0 then
          r = 128 + (r - 128) * con
          g = 128 + (g - 128) * con
          b = 128 + (b - 128) * con
          if r < 0 then r = 0 elseif r > 255 then r = 255 end
          if g < 0 then g = 0 elseif g > 255 then g = 255 end
          if b < 0 then b = 0 elseif b > 255 then b = 255 end
        end

        -- fade (0=black, 1=normal, 2=white)
        if fade < 1.0 then
          r = r * fade; g = g * fade; b = b * fade
        elseif fade > 1.0 then
          local ft = fade - 1.0
          r = r + (255 - r) * ft; g = g + (255 - g) * ft; b = b + (255 - b) * ft
        end

        pal(i, (flr(r) << 16) | (flr(g) << 8) | flr(b), 2)
      end
    end
    this._pal_dirty = true
  elseif this._pal_dirty then
    for i = 0, 63 do
      if base_pal[i] then
        pal(i, base_pal[i], 2)
      end
    end
    this._pal_dirty = false
  end
end

metaplanes.draw = function(this)
  local t = this.t

  -- draw rotating background
  if bg_tex then
    palt(0, false)
    poke2(0x5530, 0)
    poke2(0x5532, 0)
    poke2(0x5534, bg_tw)
    poke2(0x5536, bg_th)

    local W, H = CFG.w, CFG.h
    local rot = t * CFG.bg_rot_speed
    local ca, sa = cos(rot), -sin(rot)
    local ca_abs, sa_abs = abs(ca), abs(sa)
    local ihw, ihh = bg_iw * 0.5, bg_ih * 0.5
    local bw = (W * 0.5) * ca_abs + (H * 0.5) * sa_abs
    local bh = (W * 0.5) * sa_abs + (H * 0.5) * ca_abs
    local zoom = max(1.0, max(bw / ihw, bh / ihh))
    local inv_z = 1.0 / zoom

    local sx_l = (-W * 0.5) * inv_z
    local sx_r = ( W * 0.5) * inv_z
    local lu_base = ihw + sx_l * ca
    local lv_base = ihh + sx_l * sa
    local ru_base = ihw + sx_r * ca
    local rv_base = ihh + sx_r * sa
    local sy0 = (-H * 0.5) * inv_z
    local lu = lu_base - sy0 * sa
    local lv = lv_base + sy0 * ca
    local ru = ru_base - sy0 * sa
    local rv = rv_base + sy0 * ca
    local du = -sa * inv_z
    local dv =  ca * inv_z
    for y = 0, H - 1 do
      _tline3d(bg_tex, 0, y, W - 1, y, lu, lv, ru, rv, 1, 1, 0)
      lu = lu + du
      lv = lv + dv
      ru = ru + du
      rv = rv + dv
    end
  else
    cls(CFG.bg_col)
  end

  -- apply rocket params to CFG
  CFG.focal = this.focal
  CFG.cam_roll_speed = this.cam_roll
  CFG.vox_scale = this.vox_scale

  -- shake: offset drawing origin
  local shake_ox, shake_oy = 0, 0
  if this.shake > 0 then
    local a = this.shake
    shake_ox = flr((irnd(2000) / 1000 - 1) * a)
    shake_oy = flr((irnd(2000) / 1000 - 1) * a)
    camera(-shake_ox, -shake_oy)
  end

  build_field(t)
  compute_and_draw_voxels_rocket(this)

  -- reset camera offset
  if shake_ox ~= 0 or shake_oy ~= 0 then
    camera(0, 0)
  end

  -- glitch post-processing: horizontal slice displacement
  if this.glitch > 0 then
    local disp = get_display()
    local amp = this.glitch * 18
    local slices = 24
    local iw, ih = CFG.w, CFG.h
    local snap = userdata("u8", iw, ih)
    snap:copy(disp, true, 0, 0, iw * ih, 0, 0, 1)
    for s = 0, slices - 1 do
      local sy0 = flr(s * ih / slices)
      local sy1 = flr((s + 1) * ih / slices)
      local sh = sy1 - sy0
      if sh > 0 then
        local off = flr((irnd(2000) / 1000 - 1) * amp)
        if irnd(9) == 0 then off += (irnd(2) * 2 - 1) * flr(amp * 0.9) end
        if off != 0 then
          for y = sy0, sy1 - 1 do
            local row_off = y * iw
            if off > 0 then
              disp:copy(snap, true, row_off, row_off + off, iw - off, 0, 0, 1)
              disp:copy(0, true, 0, row_off, off, 0, 0, 1)
            else
              local a = -off
              disp:copy(snap, true, row_off + a, row_off, iw - a, 0, 0, 1)
              disp:copy(0, true, 0, row_off + iw - a, a, 0, 0, 1)
            end
          end
        end
      end
    end
  end
end