// Retarded-time sum over the plate for one output sample (Rayleigh integral).
// lw[k] = (weight, delay in samples) for node k; delay is relative to the
// closest node so that the ring buffer only has to cover the spread.
// The ring buffer holds the transverse velocity at the sample instants; the
// acceleration is the first difference of the delayed velocity, which is the
// average acceleration over one sample interval (band-limited, and immune to
// the per-step rounding noise of the force evaluation).

@group(0) @binding(0) var<uniform> P: Params;
@group(0) @binding(1) var<uniform> SP: StepParams;
@group(0) @binding(2) var<storage, read> hist: array<f32>;
@group(0) @binding(3) var<storage, read> lw: array<vec2f>;
@group(0) @binding(4) var<storage, read_write> audioOut: array<f32>;

var<workgroup> sh: array<f32, 256>;

@compute @workgroup_size(256)
fn main(@builtin(local_invocation_id) lid: vec3u) {
  let t = lid.x;
  let H = i32(P.histLen);
  let slot = SP.audioSlot;
  var acc = 0.0;
  for (var k = t; k < P.n; k += 256u) {
    let d = lw[k].y;
    let k0 = i32(floor(d));
    let f = d - f32(k0);
    let s0 = ((slot - k0) % H + H) % H;
    let s1 = ((slot - k0 - 1) % H + H) % H;
    let s2 = ((slot - k0 - 2) % H + H) % H;
    let v0 = hist[u32(s0) * P.n + k];
    let v1 = hist[u32(s1) * P.n + k];
    let v2 = hist[u32(s2) * P.n + k];
    let a = (mix(v0, v1, f) - mix(v1, v2, f)) * P.fs;
    acc += lw[k].x * a;
  }
  sh[t] = acc;
  workgroupBarrier();
  for (var stride = 128u; stride > 0u; stride >>= 1u) {
    if (t < stride) { sh[t] += sh[t + stride]; }
    workgroupBarrier();
  }
  if (t == 0u) { audioOut[u32(SP.audioOut)] = sh[0]; }
}