-------------------------
-- USER PARAMS
-------------------------

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

p={
  vid=3, bg=0,

  -- camera / projection
  camz=220,
  fov=210,
  near=10,
  min_s=0.14,

  -- base helix geometry (scaled per helix)
  r=64,
  pitch=7.0,
  beads=64,
  step=1,
  twist_k=0.024,     -- u -> angle slope (controls how tight helix twists along length)

  -- base animation
  time_scale=1.0,
  rots=0.0035,      -- twist speed
  cache_period=1,   -- recompute nodes every N frames (1 = every frame)
  flow_speed=0.12,  -- node flow along helix (units per frame)

  -- rung config
  rods=1,
  rod_every=3,

  -- sphere sizes in world units (scaled per helix, then projected)
  node_r=12.8,
  chain_r=3.7,      -- bead-chain radius (world)

  -- bead-chain packing (screen space)
  chain_max=16,     -- clamp chain beads per segment (stabilizes CPU)
  chain_step_scale=0.8, -- >1 = fewer beads between nodes, <1 = more beads
  rod_chain_step_scale=0.1, -- rod bead spacing scale (smaller = tighter)
  chain_jitter=3.35,     -- organic scatter for in-between beads (0..1)

  -- z buckets (used for depth shading only)
  zbins=8,

  -- sprite image
  spr_path="sphere.pod",
  spr_pal_path="sphere_00.pod.pal",
  spr_pal2_path="sphere_01.pod.pal",
  pal_fade_speed=0.3, -- full cycles per second (0.3 = ~3.3s per transition)
  spr_preload_pal=true,
  spr_apply_pal=true,
  spr_size=64,      -- source sprite size (pixels)
  spr_cols=2,       -- sprite variants per row (2x2 grid in 128x128)
  spr_rows=2,
  shade_near=7,     -- brightness shift at near bins (positive = brighter)
  shade_far=-30,     -- brightness shift at far bins (negative = darker)
  profile=false,    -- show per-section CPU deltas
  udtest=false,    -- show userdata offset test
  filter_vis=false, -- cull behind/too-small beads before sort
  use_matmul3d=true, -- use matmul3d + manual translation
  sort_mode="bin", -- "bin" (fast) or "full" (accurate)
  sort_bins=64,     -- bins for sort_mode="bin"

  -- 3 helixes
  -- ox,oy   : screen offset in pixels
  -- zoff    : world z offset (depth)
  -- sc      : uniform scale
  -- phase   : random phase (also de-correlates wobble)
  -- rot0    : initial TWIST offset
  -- rot_spd : twist speed multiplier
  -- pitch0/yaw0/roll0 : base orientation (THIS fixes “tilted away”)
  -- yaw_wob_amp/yaw_wob_spd : optional sway
  hx={
    {ox=-72, oy=  2, zoff= 70, sc=0.82, phase=1.30, rot0=0.80, rot_spd=1.00,
     pitch0=0.05, yaw0= 0.10, roll0=0.02, yaw_wob_amp=1.0, yaw_wob_spd=1.0},

    {ox=  0, oy=  0, zoff=  0, sc=1.00, phase=0.00, rot0=0.00, rot_spd=1.00,
     pitch0=0.09, yaw0= 0.00, roll0=0.04, yaw_wob_amp=1.5, yaw_wob_spd=2.0},

    {ox= 76, oy= -2, zoff= 130, sc=0.74, phase=2.10, rot0=1.60, rot_spd=0.92,
     pitch0=-0.03, yaw0=-0.08, roll0=0.01, yaw_wob_amp=1.1, yaw_wob_spd=0.85},
  },
}

dnahelix = {}

-- rocket-controlled state
local _fade_val = 1.0
local _shake_amt = 0
local _glitch_amt = 0
local _pal_blend = -1  -- -1 = auto ping-pong
local _pal_dirty = false

-- per-helix override tables (plain Lua, avoids userdata indexing issues)
local _hx_live = { {}, {}, {} }

-------------------------
-- rocket param layout (params[2..16]) — matches .kprj track order
-------------------------
-- [2]  time
-- [3]  rot_speed (global multiplier, scales each helix's own rot_spd)
-- [4]  H1 X   [5]  H2 X   [6]  H3 X
-- [7]  H1 Pitch  [8]  H2 Pitch  [9]  H3 Pitch
-- [10] H1 Yaw  [11] H2 Yaw  [12] H3 Yaw
-- [13] fade (0=black 1=normal 2=white)
-- [14] shake
-- [15] glitch
-- [16] pal_blend (0=palette A, 1=palette B, fractional for blend)

function rot_x(y,z,a) local ca,sa=cos(a),sin(a) return y*ca-z*sa, y*sa+z*ca end
function rot_y(x,z,a) local ca,sa=cos(a),sin(a) return x*ca-z*sa, x*sa+z*ca end
function rot_z(x,y,a) local ca,sa=cos(a),sin(a) return x*ca-y*sa, x*sa+y*ca end

int_cfg_ud=nil
INT_BEADS=0
INT_ROD_EVERY=1
INT_CHAIN_MAX=2
INT_SPR_SIZE=3

function init_int_cfg_ud()
  int_cfg_ud=userdata("u16", 1, 4)
  int_cfg_ud:set(0,INT_BEADS, p.beads)
  int_cfg_ud:set(0,INT_ROD_EVERY, p.rod_every)
  int_cfg_ud:set(0,INT_CHAIN_MAX, p.chain_max)
  int_cfg_ud:set(0,INT_SPR_SIZE, p.spr_size)
end

function int_get(i)
  return int_cfg_ud:get(0,i)
end

udtest_msg=nil
ud_offsets_ok=false
ud_row_major=nil

function ud_selftest()
  local w,h=5,3
  local zero=userdata("f64", w, h)
  local a=userdata("f64", w, h)
  a:set(2,1, 123)

  local found_off1=nil
  for off=0,w*h-1 do
    local b=userdata("f64", w, h)
    zero:add(a, b, off, 0, 1)
    if b:get(0,0)==123 then
      found_off1=off
      break
    end
  end

  if not found_off1 then
    ud_offsets_ok=false
    ud_row_major=nil
    udtest_msg="ud offsets unsupported"
    return
  end

  if found_off1 == (1*w + 2) then ud_row_major=true
  elseif found_off1 == (2*h + 1) then ud_row_major=false
  else ud_row_major=nil end

  local found_off2=nil
  local target_x,target_y=3,2
  for off=0,w*h-1 do
    local b=userdata("f64", w, h)
    zero:add(a, b, found_off1, off, 1)
    if b:get(target_x,target_y)==123 then
      found_off2=off
      break
    end
  end

  ud_offsets_ok = (found_off2~=nil)
  if not ud_offsets_ok then
    udtest_msg="ud offsets fail (offset2)"
  else
    local layout = ud_row_major==true and "row" or (ud_row_major==false and "col" or "unk")
    udtest_msg="ud offsets ok ("..layout..")"
  end
end

function rand01(seed)
  seed = (seed*1103515245 + 12345) & 0x7fffffff
  return seed, (seed & 0xffff) / 0x10000
end

function set_node(buf,i,x,y,z)
  local row=i-1
  buf:set(0,row,x)
  buf:set(1,row,y)
  buf:set(2,row,z)
end

function get_node(buf,i)
  local row=i-1
  return buf:get(0,row), buf:get(1,row), buf:get(2,row)
end

function get_helix_params(hi)
  -- try Lua table first (populated by setparams each frame)
  local h = _hx_live[hi]
  if h and h.sc then
    return h.ox, h.oy, h.zoff, h.sc, h.phase,
           h.rot0, h.rot_spd, h.pitch, h.yaw, h.roll,
           h.yaw_wob_amp, h.yaw_wob_spd
  end

  -- fallback: read from userdata (before first setparams call)
  local idx=hi
  local ox=helix_params:get(1,idx)
  local oy=helix_params:get(2,idx)
  local zoff=helix_params:get(3,idx)
  local sca=helix_params:get(4,idx)
  local ph=helix_params:get(5,idx)

  local rot0=helix_params:get(6,idx)
  local rot_spd=helix_params:get(7,idx)
  local pitch_a=helix_params:get(8,idx)
  local yaw_a_base=helix_params:get(9,idx)
  local roll_a=helix_params:get(10,idx)
  local yaw_wob_amp=helix_params:get(11,idx)
  local yaw_wob_spd=helix_params:get(12,idx)

  if sca==0 then
    local hp=p.hx[hi]
    if hp then
      ox=hp.ox
      oy=hp.oy
      zoff=hp.zoff
      sca=hp.sc
      ph=hp.phase
      rot0=hp.rot0
      rot_spd=hp.rot_spd
      pitch_a=hp.pitch0
      yaw_a_base=hp.yaw0
      roll_a=hp.roll0
      yaw_wob_amp=hp.yaw_wob_amp
      yaw_wob_spd=hp.yaw_wob_spd
    end
  end

  return ox or 0, oy or 0, zoff or 0, sca or 1, ph or 0,
         rot0 or 0, rot_spd or 1, pitch_a or 0, yaw_a_base or 0, roll_a or 0,
         yaw_wob_amp or 0, yaw_wob_spd or 0
end

function chain_seed(seed_id)
  return (seed_id*1103515245 + 12345) & 0x7fffffff
end

prof_labels={"build","_","_","_","xform","sort","draw"}
prof_acc={}
prof_t={}
prof_on=false

function prof_reset()
  if not prof_on then return end
  for i=1,#prof_labels do prof_acc[i]=0 end
end

function prof_begin(i)
  if not prof_on then return end
  prof_t[i]=stat(1)
end

function prof_end(i)
  if not prof_on then return end
  prof_acc[i]=(prof_acc[i] or 0) + (stat(1)-prof_t[i])
end

function prof_draw(y)
  if not prof_on then return y end
  for i=1,#prof_labels do
    local v=prof_acc[i] or 0
    v=flr(v*1000)/1000
    print(prof_labels[i]..":"..tostr(v), 2, y, 6)
    y+=7
  end
  return y
end

-- map z to shading bin using fixed range
function zbin(z)
  local b=flr((z - zmin)*zscale)+1
  if b<1 then return 1 end
  if b>zbins then return zbins end
  return b
end

-------------------------
-- palette + sprite load
-------------------------
function pack_rgba(r,g,b,a)
  a=a or 255
  return (a<<24) | (r<<16) | (g<<8) | b
end

function parse_hex_palette(text)
  if type(text)~="string" then return nil end
  local pal={}
  for token in text: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
    token=token:gsub("#","")
    local v=tonumber(token,16)
    if v then
      local r=(v>>16)&0xff
      local g=(v>>8)&0xff
      local b=v&0xff
      pal[#pal+1]={r,g,b}
    end
  end
  if #pal==0 then
    for token in text: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][0-9a-fA-F][0-9a-fA-F]") do
      token=token:gsub("#","")
      token=token:sub(-6)
      local v=tonumber(token,16)
      if v then
        local r=(v>>16)&0xff
        local g=(v>>8)&0xff
        local b=v&0xff
        pal[#pal+1]={r,g,b}
      end
    end
  end
  if #pal==0 then return nil end
  return pal
end

function load_palette(path)
  if not path or path=="" then return nil end
  local obj=fetch(path)
  -- pod.pal files: text with 0xRRGGBB values
  if type(obj)=="table" then
    obj=table.concat(obj,"\n")
  end
  if type(obj)=="string" then
    -- try 0xRRGGBB format first (pod.pal)
    local colors={}
    for token in obj:gmatch("0x([0-9a-fA-F]+)") do
      local v=tonumber(token,16)
      if v then
        local r=(v>>16)&0xff
        local g=(v>>8)&0xff
        local b=v&0xff
        colors[#colors+1]={r,g,b}
      end
    end
    if #colors>0 then return colors end
    -- fallback: #RRGGBB / RRGGBB hex format
    return parse_hex_palette(obj)
  end
  return nil
end

function apply_palette_table(ptbl)
  if not ptbl then return end
  for i=0,63 do
    local v=ptbl[i+1]
    if v then
      local r=v.r or v[1] or 0
      local g=v.g or v[2] or 0
      local b=v.b or v[3] or 0
      local a=v.a or v[4] or 255
      local packed=pack_rgba(r,g,b,a)
      pal(i,packed,2)
--      for bank=0,3 do
--        poke4(0x5000+bank*0x100+i*4, packed)
--      end
    end
  end
--  memset(0x5400,0,120)
end

-- palette fade state
pal_a=nil   -- RGB table for palette A
pal_b=nil   -- RGB table for palette B
pal_cur=nil -- current blended RGB table
pal_fade_t=0
pal_last_bin_t=-1  -- last fade value used for shade cache rebuild

function lerp_palettes(pa, pb, t)
  local out={}
  local t1=1-t
  for i=1,64 do
    local a=pa[i]
    local b=pb[i]
    local ar,ag,ab = a[1],a[2],a[3]
    local br,bg,bb = b[1],b[2],b[3]
    out[i]={flr(ar*t1+br*t+0.5), flr(ag*t1+bg*t+0.5), flr(ab*t1+bb*t+0.5)}
  end
  return out
end

function apply_and_rebuild(blend)
  pal_cur=blend
  apply_palette_table(blend)
  build_luma_ranks(blend)
  -- rebuild shade cache
  shade_cache={}
  for b=1,zbins do
    local shift=depth_shift(b)
    local remap={}
    for i=1,63 do
      local r=rank_ud:get(0,i)+shift
      if r<=0 then
        remap[i]=order_ud:get(0,0)
      elseif r>=63 then
        remap[i]=order_ud:get(0,63)
      else
        local r0=flr(r)
        local f=r-r0
        local c0=order_ud:get(0,r0)
        local c1=order_ud:get(0,r0+1)
        if ((i+b)%4)<flr(f*4) then remap[i]=c1 else remap[i]=c0 end
      end
    end
    shade_cache[b]=remap
  end
end

function is_u8_image(obj)
  if not obj then return false end
  local ok = pcall(function() return obj:width() end)
  return ok
end

function find_bmp(obj)
  if is_u8_image(obj) then return obj end
  if type(obj)~="table" then return nil end
  if is_u8_image(obj.bmp) then return obj.bmp end
  if is_u8_image(obj.image) then return obj.image end
  if is_u8_image(obj.img) then return obj.img end
  if type(obj[1])=="table" and is_u8_image(obj[1].bmp) then return obj[1].bmp end
  return nil
end

lum_ud=nil
order_ud=nil
rank_ud=nil

function ensure_luma_tables()
  if not lum_ud then lum_ud=userdata("u8", 1, 64) end
  if not order_ud then order_ud=userdata("u8", 1, 64) end
  if not rank_ud then rank_ud=userdata("u8", 1, 64) end
end

function build_luma_ranks(pal)
  ensure_luma_tables()
  for i=0,63 do
    local v=pal[i+1]
    local r=v.r or v[1] or 0
    local g=v.g or v[2] or 0
    local b=v.b or v[3] or 0
    lum_ud:set(0,i, (r*30 + g*59 + b*11)\100)
    order_ud:set(0,i, i)
  end

  for i=0,63 do
    local best=i
    local bestv=lum_ud:get(0, order_ud:get(0,i))
    for j=i+1,63 do
      local v=lum_ud:get(0, order_ud:get(0,j))
      if v<bestv then best=j bestv=v end
    end
    if best!=i then
      local a=order_ud:get(0,i)
      local b=order_ud:get(0,best)
      order_ud:set(0,i, b)
      order_ud:set(0,best, a)
    end
  end

  for i=0,63 do
    rank_ud:set(0, order_ud:get(0,i), i)
  end
end

-------------------------
-------------------------
-- palette remap by depth
-------------------------
function depth_shift(bin)
  local denom = max(1, (zbins-1))
  local t = (zbins - bin) / denom
  -- smoothstep for a softer ramp
  t = t*t*(3-2*t)
  return flr(p.shade_far + (p.shade_near - p.shade_far) * t)
end

function set_depth_shade_cached(bin)
  pal(0)
  pal(0,0)
  palt(0,true)
  local c=shade_cache[bin]
  if not c then return end
  for i=1,63 do pal(i,c[i]) end
end

function set_depth_shade(bin, shift)
  pal(0)
  if shift==nil then shift=depth_shift(bin) end
  pal(0,0) -- keep color 0 transparent
  palt(0,true)
  for i=1,63 do
    local r=rank_ud:get(0,i) + shift
    if r<=0 then
      pal(i, order_ud:get(0,0))
    elseif r>=63 then
      pal(i, order_ud:get(0,63))
    else
      local r0=flr(r)
      local f=r - r0
      local c0=order_ud:get(0,r0)
      local c1=order_ud:get(0,r0+1)
      -- ordered dithering between neighboring ranks for smoother shading
      if ((i + bin) % 4) < flr(f*4) then
        pal(i, c1)
      else
        pal(i, c0)
      end
    end
  end
end

-------------------------
-- prealloc / state
-------------------------
function dnahelix.load(_self)
  if not int_cfg_ud then init_int_cfg_ud() end
  if vid then vid(p.vid) end
  w,h=480,270
  if p.vid==3 then w,h=240,135 end
  if p.vid==4 then w,h=160,90 end
  cx,cy=w*0.5,h*0.5

  -- load palette + sprite image
  local pal1=nil
  local pal2=nil
  if p.spr_apply_pal then
    pal1=load_palette(p.spr_pal_path)
    pal2=load_palette(p.spr_pal2_path)
    if p.spr_preload_pal and pal1 then
      apply_palette_table(pal1)
    end
  end

  local obj=fetch(p.spr_path)
  sprite=unpod(obj)
  if not sprite then sprite=find_bmp(obj) end
  if not sprite then
    sprite=userdata("u8", int_get(INT_SPR_SIZE), int_get(INT_SPR_SIZE))
  end

  if p.spr_apply_pal and pal1 then
    -- store both palettes for fading
    pal_a=pal1
    pal_b=pal2 or pal1  -- fallback to same if pal2 missing
    pal_cur=pal1
    apply_palette_table(pal1)
    build_luma_ranks(pal1)
  else
    -- fallback: build ranks from current palette
    local cur={}
    for i=0,63 do
      local v=peek4(0x5000+i*4)
      local r=(v>>16)&0xff
      local g=(v>>8)&0xff
      local b=v&0xff
      cur[i+1]={r,g,b}
    end
    build_luma_ranks(cur)
  end

  -- reusable node buffers (userdata)
  strand_a_nodes=userdata("f32", 3, int_get(INT_BEADS))
  strand_b_nodes=userdata("f32", 3, int_get(INT_BEADS))

  -- estimate max sprite draws per helix:
  -- nodes: 2*beads
  -- strand segments: 2*(beads-1) each with <=chain_max beads
  -- rungs: ~beads/rod_every each with <=chain_max beads
  local nh=#p.hx
  local n=int_get(INT_BEADS)
  local chain_max=int_get(INT_CHAIN_MAX)
  local rung_ct = (p.rods==1) and ceil(n/int_get(INT_ROD_EVERY)) or 0
  max_beads_per_helix = 2*n + 2*(n-1)*chain_max + rung_ct*chain_max
  max_beads = max_beads_per_helix * nh + 16

  -- helix params in userdata (per-helix)
  -- cols: 1=ox,2=oy,3=zoff,4=sc,5=phase,6=rot0,7=rot_spd,8=pitch0,9=yaw0,10=roll0,11=yaw_wob_amp,12=yaw_wob_spd
  -- NOTE: uses 1-based row indexing (row 0 unused), so need nh+1 rows; col 12 needs width 13
  helix_params=userdata("f32", 13, nh+1)
  for i=1,nh do
    local h0=p.hx[i]
    local r=i
    helix_params:set(1,r, h0.ox)
    helix_params:set(2,r, h0.oy)
    helix_params:set(3,r, h0.zoff)
    helix_params:set(4,r, h0.sc)
    helix_params:set(5,r, h0.phase)
    helix_params:set(6,r, h0.rot0)
    helix_params:set(7,r, h0.rot_spd)
    helix_params:set(8,r, h0.pitch0)
    helix_params:set(9,r, h0.yaw0)
    helix_params:set(10,r, h0.roll0)
    helix_params:set(11,r, h0.yaw_wob_amp)
    helix_params:set(12,r, h0.yaw_wob_spd)
  end

  zbins=p.zbins

  -- center pixel color for tiny-bead rendering (circfill/pset)
  local half_spr=int_get(INT_SPR_SIZE)\2
  spr_mid_col=sprite:get(half_spr, half_spr)

  -- pre-compute shade palette remaps per depth bin
  shade_cache={}
  for b=1,zbins do
    local shift=depth_shift(b)
    local remap={}
    for i=1,63 do
      local r=rank_ud:get(0,i)+shift
      if r<=0 then
        remap[i]=order_ud:get(0,0)
      elseif r>=63 then
        remap[i]=order_ud:get(0,63)
      else
        local r0=flr(r)
        local f=r-r0
        local c0=order_ud:get(0,r0)
        local c1=order_ud:get(0,r0+1)
        if ((i+b)%4)<flr(f*4) then remap[i]=c1 else remap[i]=c0 end
      end
    end
    shade_cache[b]=remap
  end

  -- helix bead buffers (world-space) stored in userdata
  -- cols: x,y,z,w(=1),size
  helix_cols=5
  helix_local=userdata("f64", helix_cols, max_beads_per_helix)
  transformed_helix=userdata("f64", helix_cols, max_beads_per_helix)
  all_helix=userdata("f64", helix_cols, max_beads)
  zero_helix=userdata("f64", helix_cols, max_beads_per_helix)
  zero_all=userdata("f64", helix_cols, max_beads)
  all_bin_idx=userdata("u32", max_beads, p.sort_bins)
  all_bin_cnt=userdata("u32", 1, p.sort_bins)
  has_matmul3d = (type(helix_local.matmul3d)=="function")
  helix_bin_idx=userdata("u32", max_beads_per_helix, zbins)
  helix_bin_cnt=userdata("u32", 1, zbins)
  helix_matrix=userdata("f64", 4, 4)
  helix_matrix_t=userdata("f64", 4, 4)
  helix_count=0

  -- z-range (cover all helixes)
  local minzoff,maxzoff=0,0
  local maxr=0
  for i=1,nh do
    local h0=p.hx[i]
    if h0.zoff<minzoff then minzoff=h0.zoff end
    if h0.zoff>maxzoff then maxzoff=h0.zoff end
    local rr=p.r*h0.sc
    if rr>maxr then maxr=rr end
  end
  zmin = (-maxr - 24) + minzoff
  zmax = ( maxr + 24) + maxzoff
  zscale = zbins / (zmax - zmin)

  ud_selftest()

  t=0
end

function dnahelix.init(_self)
  t=0
  pal_fade_t=0
  pal_last_bin_t=-1
  _fade_val=1.0
  _shake_amt=0
  _glitch_amt=0
  _pal_blend=-1
  _pal_dirty=false
end

function dnahelix.setparams(_self, params)
  -- time (rocket controls animation time directly)
  if params[2] ~= nil then t = params[2] end

  -- global rotation speed multiplier (scales each helix's own rot_spd)
  -- rocket sends 0 for unset tracks, so default to 1.0 when 0
  local global_rot = params[3]
  if global_rot == nil or global_rot == 0 then global_rot = 1.0 end

  -- per-helix params: project layout groups by property, not by helix
  -- params[4,5,6] = H1/H2/H3 X
  -- params[7,8,9] = H1/H2/H3 Pitch
  -- params[10,11,12] = H1/H2/H3 Yaw
  for hi=1,3 do
    local h0 = p.hx[hi]
    if not h0 then goto next_helix end

    local r_ox     = h0.ox     + (params[3 + hi] or 0)   -- params[4,5,6]
    local r_pitch  = h0.pitch0 + (params[6 + hi] or 0) * 0.01   -- params[7,8,9] scaled: 25=quarter turn
    local r_yaw    = h0.yaw0   + (params[9 + hi] or 0) * 0.01   -- params[10,11,12] scaled: 25=quarter turn
    local r_rot_spd = h0.rot_spd * global_rot

    -- store in plain Lua table (reliable)
    _hx_live[hi] = {
      ox   = r_ox,
      oy   = h0.oy,
      zoff = h0.zoff,
      sc   = h0.sc,
      phase = h0.phase,
      rot0 = h0.rot0,
      rot_spd = r_rot_spd,
      pitch = r_pitch,
      yaw  = r_yaw,
      roll = h0.roll0,
      yaw_wob_amp = h0.yaw_wob_amp,
      yaw_wob_spd = h0.yaw_wob_spd,
    }

    -- also write to userdata (legacy)
    helix_params:set(1,  hi, r_ox)
    helix_params:set(8,  hi, r_pitch)
    helix_params:set(9,  hi, r_yaw)
    helix_params:set(2,  hi, h0.oy)
    helix_params:set(3,  hi, h0.zoff)
    helix_params:set(4,  hi, h0.sc)
    helix_params:set(5,  hi, h0.phase)
    helix_params:set(6,  hi, h0.rot0)
    helix_params:set(7,  hi, r_rot_spd)
    helix_params:set(10, hi, h0.roll0)
    helix_params:set(11, hi, h0.yaw_wob_amp)
    helix_params:set(12, hi, h0.yaw_wob_spd)
    ::next_helix::
  end

  -- global controls (0 = safe default for all three)
  _fade_val   = params[13] == nil and 1.0 or params[13]
  _shake_amt  = params[14] or 0
  _glitch_amt = params[15] or 0
  _pal_blend  = params[16] == nil and -1 or params[16]
end

function dnahelix.update(_self)
  -- advance time each frame (rocket overrides via setparams)
  t += 1

  -- palette blend (rocket-controlled or auto ping-pong)
  if pal_a and pal_b then
    local blend_t
    if _pal_blend >= 0 then
      -- rocket-controlled: 0=pal_a, 1=pal_b
      blend_t = _pal_blend
      if blend_t > 1 then blend_t = 1 end
    else
      -- auto ping-pong fallback
      pal_fade_t += p.pal_fade_speed / 60
      local phase = pal_fade_t % 2
      blend_t = phase <= 1 and phase or (2 - phase)
      blend_t = blend_t * blend_t * (3 - 2 * blend_t)
    end
    -- quantize to ~32 steps to avoid rebuilding shade cache every frame
    local qt = flr(blend_t * 32 + 0.5) / 32
    if qt ~= pal_last_bin_t then
      pal_last_bin_t = qt
      local blended = lerp_palettes(pal_a, pal_b, qt)
      apply_and_rebuild(blended)
    end
  end
end

-------------------------
-- emit bead into userdata list
-------------------------
function emit_bead(buf, count, max_count, x,y,z, rr, seed_id, fixed_variant)
  if rr<0.5 then return count end
  if count>=max_count then return count end

  local row=count
  buf:set(0,row, x)
  buf:set(1,row, y)
  buf:set(2,row, z)
  buf:set(3,row, 1)
  buf:set(4,row, rr)
  return count+1
end

-- emit bead chain between two world-space endpoints
-- z is interpolated for sorting
function emit_chain(buf, count, max_count, x0,y0,z0, x1,y1,z1, rr_mid, seed_id, jitter_scale, step_scale, fixed_variant,
                    _camz, _near, _fov, _chain_step_scale, _chain_max, _chain_jitter, _cx, _cy, _w, _h)
  local dx=x1-x0
  local dy=y1-y0
  local dz=z1-z0
  local L2=dx*dx+dy*dy+dz*dz
  if L2<=0 then return count end
  -- use locals passed from caller (fall back to globals for compat)
  local _cz = _camz or camz
  local _nr = _near or near
  local _fv = _fov or fov
  local _css = _chain_step_scale or chain_step_scale
  local _cm = _chain_max or chain_max
  local _cj = _chain_jitter or chain_jitter
  -- early depth cull (before any heavy work)
  local dz_mid=(z0+z1)*0.5 + _cz
  if dz_mid<_nr then return count end
  local s_mid=_fv/dz_mid
  local L=L2^0.5
  -- LOD: skip tiny screen-space segments, cap bead count
  local seg_px=L*s_mid
  if seg_px<2 then return count end
  -- lateral frustum cull (projected midpoint +/- projected half-length)
  if _cx then
    local mx=(x0+x1)*0.5*s_mid+_cx
    local my=(y0+y1)*0.5*s_mid+_cy
    local half_px=seg_px*0.5+rr_mid*s_mid
    if mx+half_px<0 or mx-half_px>_w or my+half_px<0 or my-half_px>_h then
      return count
    end
  end
  local step=max(2, rr_mid*2*_css*(step_scale or 1))
  local ceil=math.ceil
  local m=ceil(L/step)
  if m<1 then m=1 end
  if m>_cm then m=_cm end
  local m_lod=max(1,ceil(seg_px*0.4))
  if m_lod<m then m=m_lod end
  local inv=1/m
  local nx,ny=0,0
  local Lxy2=dx*dx+dy*dy
  if Lxy2>0 then
    local invL=1/Lxy2^0.5
    nx=-dy*invL
    ny= dx*invL
  end
  local base_seed=chain_seed(seed_id or 0)
  local jitter=0
  if s_mid>0 then
    jitter=(_cj*rr_mid*(jitter_scale or 1)) / s_mid
  end
  if jitter>0 then
    -- jittered path: single hash -> 3 random values via bit extraction
    local cj015=_cj*0.15
    local rr_lo=rr_mid*0.85
    local rr_sc=rr_mid*0.3
    local inv1024=0.0009765625
    local buf_set=buf.set
    for i=1,m-1 do
      if count>=max_count then break end
      local u=i*inv
      local h=(base_seed^^(i*374761393))
      h=(h*1103515245+12345)&0x7fffffff
      local r1=(h&0x3ff)*inv1024
      local r2=((h>>10)&0x3ff)*inv1024
      local r3=((h>>20)&0x3ff)*inv1024
      u += (r1-0.5)*cj015
      if u<0 then u=0 elseif u>1 then u=1 end
      local off=(r2-0.5)*2*jitter
      local rr=rr_lo+rr_sc*r3
      if rr>=0.5 then
        buf_set(buf,0,count, x0+dx*u+nx*off)
        buf_set(buf,1,count, y0+dy*u+ny*off)
        buf_set(buf,2,count, z0+dz*u)
        buf_set(buf,3,count, 1)
        buf_set(buf,4,count, rr)
        count+=1
      end
    end
  else
    -- no-jitter path (rung chains): skip rand entirely
    local buf_set=buf.set
    for i=1,m-1 do
      if count>=max_count then break end
      local u=i*inv
      buf_set(buf,0,count, x0+dx*u)
      buf_set(buf,1,count, y0+dy*u)
      buf_set(buf,2,count, z0+dz*u)
      buf_set(buf,3,count, 1)
      buf_set(buf,4,count, rr_mid)
      count+=1
    end
  end
  return count
end

function swap_ud_rows(ud, cols, i, j)
  for c=0,cols-1 do
    local tmp=ud:get(c,i)
    ud:set(c,i, ud:get(c,j))
    ud:set(c,j, tmp)
  end
end

function sort_ud_by_z(ud, cols, n)
  local gap=flr(n/2)
  while gap>0 do
    for i=gap,n-1 do
      local j=i
      while j>=gap do
        local z0=ud:get(3,j-gap)
        local z1=ud:get(3,j)
        if z0>=z1 then break end
        swap_ud_rows(ud, cols, j, j-gap)
        j-=gap
      end
    end
    gap=flr(gap/2)
  end
end

function find_visible_start(ud, n, near)
  local lo=0
  local hi=n-1
  local ans=-1
  while lo<=hi do
    local mid=flr((lo+hi)/2)
    local dz=ud:get(3,mid)
    if dz>=near then
      ans=mid
      lo=mid+1
    else
      hi=mid-1
    end
  end
  return ans
end

function build_matrix(m, yaw, pitch, roll, tx, ty, tz)
  local cz,sz=cos(roll),sin(roll)
  local cx,sx=cos(pitch),sin(pitch)
  local cy,sy=cos(yaw),sin(yaw)

  local r00 = cy*cz + sy*sx*sz
  local r01 = -cy*sz + sy*sx*cz
  local r02 = sy*cx
  local r10 = cx*sz
  local r11 = cx*cz
  local r12 = -sx
  local r20 = -sy*cz + cy*sx*sz
  local r21 = sy*sz + cy*sx*cz
  local r22 = cy*cx

  m:set(0,0, r00) m:set(1,0, r01) m:set(2,0, r02) m:set(3,0, tx)
  m:set(0,1, r10) m:set(1,1, r11) m:set(2,1, r12) m:set(3,1, ty)
  m:set(0,2, r20) m:set(1,2, r21) m:set(2,2, r22) m:set(3,2, tz)
  m:set(0,3, 0)   m:set(1,3, 0)   m:set(2,3, 0)   m:set(3,3, 1)
end

function matmul_fallback(src, dst, count, m)
  for i=0,count-1 do
    local x=src:get(0,i)
    local y=src:get(1,i)
    local z=src:get(2,i)
    local w=src:get(3,i)
    local nx = m:get(0,0)*x + m:get(1,0)*y + m:get(2,0)*z + m:get(3,0)*w
    local ny = m:get(0,1)*x + m:get(1,1)*y + m:get(2,1)*z + m:get(3,1)*w
    local nz = m:get(0,2)*x + m:get(1,2)*y + m:get(2,2)*z + m:get(3,2)*w
    dst:set(0,i, nx)
    dst:set(1,i, ny)
    dst:set(2,i, nz)
    dst:set(3,i, 1)
  end
end

-------------------------
-- main draw
-------------------------
function dnahelix.draw(_self)
  cls(p.bg)

  -- shake: offset camera
  local shk = _shake_amt
  local ox, oy = 0, 0
  if shk > 0 then
    ox = irnd(flr(shk * 2 + 1)) - flr(shk)
    oy = irnd(flr(shk * 2 + 1)) - flr(shk)
    camera(-ox, -oy)
  end

  prof_on=p.profile
  prof_reset()

  local l_camz=p.camz
  local l_fov=p.fov
  local l_near=p.near
  local l_time_scale=p.time_scale
  -- keep globals for backward compat with emit_chain fallback
  camz=l_camz
  fov=l_fov
  near=l_near
  time_scale=l_time_scale

  local l_nh=#p.hx
  nh=l_nh

  local l_chain_step_scale=p.chain_step_scale
  local l_chain_max=p.chain_max
  local l_chain_jitter=p.chain_jitter
  local l_min_s=p.min_s
  local l_rods=p.rods
  local l_rod_every=int_get(INT_ROD_EVERY)
  chain_step_scale=l_chain_step_scale
  chain_max=l_chain_max
  chain_jitter=l_chain_jitter
  min_s=l_min_s
  rods=l_rods
  rod_every=l_rod_every

  local l_spr=sprite
  local l_spr_size=int_get(INT_SPR_SIZE)
  local l_sspr_fn=sspr
  spr_size=l_spr_size
  sspr_fn=l_sspr_fn

  local l_strand_a=strand_a_nodes
  local l_strand_b=strand_b_nodes
  strand_a_nodes_t=l_strand_a
  strand_b_nodes_t=l_strand_b

  local n=int_get(INT_BEADS)
  local sin_fn=sin
  local cos_fn=cos
  local emit_chain_fn=emit_chain
  local build_matrix_fn=build_matrix
  local matmul_fn=matmul_fallback
  local use_matmul3d = p.use_matmul3d and has_matmul3d
  local sort_fn=sort_ud_by_z
  local find_start_fn=find_visible_start
  local zbin_fn=zbin
  local shade_fn=shade_cache and set_depth_shade_cached or set_depth_shade
  local pset_fn=pset
  local cfill_fn=circfill
  local mid_col=spr_mid_col
  local use_udops = ud_offsets_ok and (ud_row_major==true) and not p.filter_vis
  local sort_mode=p.sort_mode
  local sort_bins=p.sort_bins
  local bin_idx=all_bin_idx
  local bin_cnt=all_bin_cnt
  local l_cx=cx
  local l_cy=cy
  local l_w=w
  local l_h=h
  local l_abs=abs
  local l_flr=flr
  local l_max_beads=max_beads_per_helix

  local l_twist_k=p.twist_k
  local l_flow_speed=p.flow_speed
  local l_rots=p.rots
  twist_k=l_twist_k
  flow_speed=l_flow_speed
  rots=l_rots
  local cache_period=p.cache_period
  local step=p.step
  local r_base=p.r
  local pitch_base=p.pitch
  local node_r_base=p.node_r
  local chain_r_base=p.chain_r
  local rod_chain_ss=p.rod_chain_step_scale

  local tt=t*l_time_scale
  local half_n = n*0.5

  helix_count=0
  local total_count=0

  -- build helixes: merged node computation + bead emission + chain emission
  for hi=1,l_nh do
    local ox,oy,zoff,sca,ph,rot0,rot_spd,pitch_a,yaw_a_base,roll_a,yaw_wob_amp,yaw_wob_spd = get_helix_params(hi)
    local radius=(r_base or p.r) * (sca or 1)
    local pitch=(pitch_base or p.pitch) * (sca or 1)
    local node_r=(node_r_base or p.node_r) * (sca or 1)
    local chain_r=(chain_r_base or p.chain_r) * (sca or 1)
    local node_r_big=node_r
    local node_r_small=node_r*(2/3)

    local a0 = tt*l_rots * (rot_spd or 1) + (rot0 or 0)
    local yaw_w = sin_fn(a0*0.21*(yaw_wob_spd or 0) + (ph or 0))*0.04*(yaw_wob_amp or 0)
    local yaw_a   = yaw_a_base + yaw_w

    local local_count=0
    local max_gap = pitch*step*1.5
    local hl=helix_local
    local hl_set=hl.set
    local hi_seed=hi*100000

    -- merged: compute nodes + emit node beads + emit strand chains + emit rungs
    -- all in a single pass with cached previous node
    prof_begin(1)
    local span=n*step
    local span_half=span*0.5
    local flow_offset=t*l_flow_speed + span_half
    local base_u=(-half_n) * step
    local prev_ax,prev_ay,prev_az
    local prev_bx,prev_by,prev_bz
    local do_rods = (l_rods==1) and (l_rod_every>0)

    for i=1,n do
      local u=((base_u + flow_offset) % span) - span_half
      local ang=a0 + ph + u*l_twist_k
      local ca=cos_fn(ang)
      local sa=sin_fn(ang)
      local ax=ca*radius
      local az=sa*radius
      local ay=u*pitch
      local bx=-ax
      local bz=-az
      local by=ay
      base_u += step

      -- store nodes for matmul (still needed for transform)
      local row=i-1
      l_strand_a:set(0,row,ax) l_strand_a:set(1,row,ay) l_strand_a:set(2,row,az)
      l_strand_b:set(0,row,bx) l_strand_b:set(1,row,by) l_strand_b:set(2,row,bz)

      -- inline emit_bead for both strands
      local has_rod = do_rods and ((i%l_rod_every)==0)
      local nr = has_rod and node_r_big or node_r_small
      if nr>=0.5 and local_count<l_max_beads then
        hl_set(hl,0,local_count,ax) hl_set(hl,1,local_count,ay) hl_set(hl,2,local_count,az)
        hl_set(hl,3,local_count,1)  hl_set(hl,4,local_count,nr)
        local_count+=1
      end
      if nr>=0.5 and local_count<l_max_beads then
        hl_set(hl,0,local_count,bx) hl_set(hl,1,local_count,by) hl_set(hl,2,local_count,bz)
        hl_set(hl,3,local_count,1)  hl_set(hl,4,local_count,nr)
        local_count+=1
      end

      -- strand chains from prev -> current (skip first node)
      if i>1 then
        local gap = l_abs(ay - prev_ay)
        if gap <= max_gap then
          local_count=emit_chain_fn(hl, local_count, l_max_beads, prev_ax,prev_ay,prev_az, ax,ay,az, chain_r, hi_seed+i, 1, 1, nil,
            l_camz, l_near, l_fov, l_chain_step_scale, l_chain_max, l_chain_jitter, l_cx, l_cy, l_w, l_h)
          local_count=emit_chain_fn(hl, local_count, l_max_beads, prev_bx,prev_by,prev_bz, bx,by,bz, chain_r, hi_seed+1000+i, 1, 1, nil,
            l_camz, l_near, l_fov, l_chain_step_scale, l_chain_max, l_chain_jitter, l_cx, l_cy, l_w, l_h)
        end
      end

      -- rung chain (A <-> B)
      if has_rod then
        local_count=emit_chain_fn(hl, local_count, l_max_beads, ax,ay,az, bx,by,bz, chain_r, hi_seed+2000+i, 0, rod_chain_ss, 0,
          l_camz, l_near, l_fov, l_chain_step_scale, l_chain_max, l_chain_jitter, l_cx, l_cy, l_w, l_h)
      end

      prev_ax,prev_ay,prev_az = ax,ay,az
      prev_bx,prev_by,prev_bz = bx,by,bz
    end
    prof_end(1)

    -- pipeline: matmul -> gather -> sort -> project -> draw
    prof_begin(5)
    build_matrix_fn(helix_matrix, yaw_a, pitch_a, roll_a, ox, oy, zoff)
    -- also store nodes back into helix_local for size column (already there from inline emit)
    if use_matmul3d then
      -- matmul3d uses row-vector convention (v*M), so transpose the
      -- column-vector matrix.  The transposed M^T naturally has
      -- (tx,ty,tz,1) in row 3, and w=1 is assumed, so translation
      -- is included for free — no separate :add calls needed.
      helix_matrix:transpose(helix_matrix_t)
      helix_local:matmul3d(helix_matrix_t, transformed_helix, 1)
    else
      matmul_fn(helix_local, transformed_helix, local_count, helix_matrix)
    end

    if use_udops then
      local w=helix_cols
      local base=total_count*w
      -- size column copy
      zero_helix:add(helix_local, transformed_helix, 4, 4, 1, w, w, local_count)
      -- bulk copy rows into all_helix
      zero_all:add(transformed_helix, all_helix, 0, base, w, w, w, local_count)
      -- dz = z + camz
      zero_all:add(all_helix, all_helix, base+2, base+3, 1, w, w, local_count)
      all_helix:add(l_camz, true, 0, base+3, 1, w, w, local_count)
    else
      -- copy into global buffer (safe per-element) with visibility cull
      local kept=0
      for i=0,local_count-1 do
        local z=transformed_helix:get(2,i)
        local dz=z + l_camz
        if dz>=l_near then
          local s=l_fov/dz
          if s>=l_min_s then
            local dst=total_count+kept
            all_helix:set(0,dst, transformed_helix:get(0,i))
            all_helix:set(1,dst, transformed_helix:get(1,i))
            all_helix:set(2,dst, z)
            all_helix:set(3,dst, dz)
            all_helix:set(4,dst, helix_local:get(4,i))
            kept+=1
          end
        end
      end
      local_count=kept
    end
    total_count += local_count
    prof_end(5)
  end

  if sort_mode=="bin" then
    -- fast bin sort (approx order)
    prof_begin(6)
    local zmin_dz = zmin + l_camz
    local zscale_sort = sort_bins / (zmax - zmin)
    for b=0,sort_bins-1 do
      bin_cnt:set(0,b,0)
    end
    for i=0,total_count-1 do
      local dz=all_helix:get(3,i)
      if dz>=l_near then
        local bin=l_flr((dz - zmin_dz)*zscale_sort)+1
        if bin<1 then bin=1 elseif bin>sort_bins then bin=sort_bins end
        local row=bin-1
        local c=bin_cnt:get(0,row)
        bin_idx:set(c,row,i)
        bin_cnt:set(0,row,c+1)
      end
    end
    prof_end(6)

    prof_begin(7)
    local ah=all_helix
    -- pre-compute shade-bin per sort-bin (avoids per-bead z read + zbin)
    local zb_ratio=zbins/sort_bins
    last_bin=-1
    for b=sort_bins-1,0,-1 do
      local c=bin_cnt:get(0,b)
      if c>0 then
        -- set shade once per sort-bin (derived from bin index)
        local sbin=l_flr((b+0.5)*zb_ratio)+1
        if sbin<1 then sbin=1 elseif sbin>zbins then sbin=zbins end
        if sbin~=last_bin then
          shade_fn(sbin)
          last_bin=sbin
        end
        for j=0,c-1 do
          local i=bin_idx:get(j,b)
          local dz=ah:get(3,i)
          local s=l_fov/dz
          if s<l_min_s then goto skip_bin_draw end
          local rr=ah:get(4,i)
          local diam=rr*s*2
          if diam<1 then goto skip_bin_draw end
          local px=ah:get(0,i)*s+l_cx
          local py=ah:get(1,i)*s+l_cy
          local half=diam*0.5
          if px+half<0 or px-half>l_w or py+half<0 or py-half>l_h then goto skip_bin_draw end
          if diam<2 then
            pset_fn(px,py,mid_col)
          elseif diam<4 then
            cfill_fn(px,py,half,mid_col)
          else
            l_sspr_fn(l_spr, 0,0, l_spr_size,l_spr_size, px-half,py-half, diam,diam, false,false)
          end
          ::skip_bin_draw::
        end
      end
    end
    prof_end(7)
  else
    -- full sort (accurate)
    prof_begin(6)
    sort_fn(all_helix, helix_cols, total_count)
    local end_i=find_start_fn(all_helix, total_count, near)
    prof_end(6)

    prof_begin(7)
    local ah=all_helix
    last_bin=-1
    for i=0,end_i do
      local dz=ah:get(3,i)
      if dz<l_near then goto skip_full_draw end
      local s=l_fov/dz
      if s<l_min_s then goto skip_full_draw end
      local rr=ah:get(4,i)
      local diam=rr*s*2
      if diam<1 then goto skip_full_draw end
      local px=ah:get(0,i)*s+l_cx
      local py=ah:get(1,i)*s+l_cy
      local half=diam*0.5
      if px+half<0 or px-half>l_w or py+half<0 or py-half>l_h then goto skip_full_draw end
      local z=ah:get(2,i)
      local bin=zbin_fn(z)
      if bin~=last_bin then
        shade_fn(bin)
        last_bin=bin
      end
      if diam<2 then
        pset_fn(px,py,mid_col)
      elseif diam<4 then
        cfill_fn(px,py,half,mid_col)
      else
        l_sspr_fn(l_spr, 0,0, l_spr_size,l_spr_size, px-half,py-half, diam,diam, false,false)
      end
      ::skip_full_draw::
    end
    prof_end(7)
  end


  pal(0)
  palt()

  -- reset shake camera
  if shk > 0 then camera(0, 0) end

  -- fade (0=black, 1=normal, 2=white) applied to display palette
  if _fade_val ~= 1.0 then
    local f = _fade_val
    for i = 0, 63 do
      -- read from stored palette table, not VRAM (avoids fight with palette blend)
      local r, g, b = 0, 0, 0
      if pal_cur then
        local v = pal_cur[i + 1]
        if v then r, g, b = v[1] or 0, v[2] or 0, v[3] or 0 end
      else
        -- fallback: read VRAM bank 0
        local c = peek4(0x5000 + i * 4)
        r = (c >> 16) & 0xff
        g = (c >> 8) & 0xff
        b = c & 0xff
      end
      if f < 1 then
        r = flr(r * f + 0.5)
        g = flr(g * f + 0.5)
        b = flr(b * f + 0.5)
      else
        local f2 = f - 1
        r = flr(r + (255 - r) * f2 + 0.5)
        g = flr(g + (255 - g) * f2 + 0.5)
        b = flr(b + (255 - b) * f2 + 0.5)
      end
      pal(i, (r << 16) | (g << 8) | b, 2)
    end
    _pal_dirty = true
  elseif _pal_dirty then
    -- restore the blended palette, not Picotron defaults
    if pal_cur then
      apply_palette_table(pal_cur)
    else
      for i = 0, 63 do pal(i, i, 2) end
    end
    _pal_dirty = false
  end

  -- glitch: horizontal slice displacement
  if _glitch_amt > 0 then
    local disp = get_display()
    local amp = _glitch_amt * 18
    local slices = 24
    local iw, ih = w, 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
