/** AudioWorklet processor: a bank of damped modal resonators driven by rigid exciters through a
 * power-law contact law, with the listener pressure as output.
 *
 * Mode m obeys  a'' + 2 sigma a' + omega0^2 a = f_m  with f_m = sum_p phi_m(p) (-F_p), the modal
 * force of the contact springs at the patch points p. Its state is a complex number z with
 * a = Im(z) and  dz/dt = lambda z + f_m / omega_d,  lambda = -sigma + i omega_d. One sample with
 * the force held constant is  z <- P z + b f_m  with P = exp(lambda dt) and
 * b = (P - 1) / (lambda omega_d), exact for piecewise-constant forcing. The real acceleration is
 * Im(lambda^2 z) + f_m, so the listener pressure with the complex per-mode gain G_m (Rayleigh
 * sum, delays folded into a phase at the mode frequency) is
 *     p = sum_m Im(G_m lambda_m^2 z_m) + Re(G_m) f_m .
 *
 * Contact is implicit: the force of a sample is the one that satisfies the contact law at the
 * end of that sample, F_p = k max(w_p + o_p - z_s, 0)^alpha with w_p and z_s advanced by the
 * force itself. The surface under the patch responds to the step's forces through the P x P
 * compliance c_pq = sum_m phi_m(p) phi_m(q) Im(b_m), and the exciter through dt^2 / m, so the
 * end-of-step penetration is affine in F and the small nonlinear system is solved by Newton
 * iteration. This keeps the bank at one step per sample for any contact stiffness, where an
 * explicit coupling would need sub-steps. Between hits the bank runs mode-outer over the block.
 *
 * A hit carries the frame at which its exciter is released, so one processor serves live play
 * (release at once), timed playback of a recording, and an offline render, where the whole
 * performance is handed over in processorOptions before rendering begins. A hit whose frame has
 * not arrived costs nothing until the block it falls in.
 */

interface InitMessage {
  type: "init";
  freq: Float64Array;
  sigma: Float64Array;
  gain: Float32Array; // (Nm, 2) complex listener gain, re/im interleaved
}

/** New frequencies and decay rates for the same modes, keeping the ringing state: this is what
 * the design panel sends while damping or material sliders move. */
interface ParamsMessage {
  type: "params";
  freq: Float64Array;
  sigma: Float64Array;
}

interface GainMessage {
  type: "gain";
  gain: Float32Array;
}

interface HitMessage {
  type: "hit";
  id: number;
  phi: Float32Array; // (P, Nm) point-major mode values at the patch points
  offsets: Float32Array; // (P,) surface height minus exciter face plane
  kPoint: number; // contact stiffness per point [N/m^alpha]
  alpha: number;
  mass: number;
  gap: number;
  velocity: number; // approach speed, positive
  atFrame?: number; // release frame; absent means as soon as possible
  silent?: boolean; // do not report this hit back (playback and offline renders)
}

type InMessage = InitMessage | ParamsMessage | GainMessage | HitMessage | { type: "reset" };

interface Exciter {
  id: number;
  phi: Float32Array;
  offsets: Float32Array;
  nPoints: number;
  kPoint: number;
  alpha: number;
  mass: number;
  gap: number;
  velocity: number;
  atFrame: number;
  silent: boolean;
  released: boolean;
  zs: number;
  vs: number;
  startFrame: number;
  peakForce: number;
  contactSamples: number;
  newtonIters: number;
  touched: boolean;
  wp: Float64Array; // free end-of-step surface height under the points (with offsets)
  Fp: Float64Array; // current forces
  c: Float64Array; // (P, P) one-step compliance of the surface
  J: Float64Array; // Newton work: (P, P)
  g: Float64Array;
  dF: Float64Array;
}

const HOLD = 1e3;
const now: () => number = typeof performance !== "undefined" ? () => performance.now() : () => Date.now();
const RELEASE_CLEARANCE = 2e-3; // an exciter is dropped once it has rebounded this far clear
const NEWTON_MAX = 12;

class CymbalModalProcessor extends AudioWorkletProcessor {
  private nModes = 0;
  private zr = new Float64Array(0);
  private zi = new Float64Array(0);
  private lamR = new Float64Array(0);
  private lamI = new Float64Array(0);
  private pr = new Float64Array(0); // one-sample pole
  private pi = new Float64Array(0);
  private br = new Float64Array(0); // one-sample forcing coefficient
  private bi = new Float64Array(0);
  private cr = new Float64Array(0); // output coefficient C = G lambda^2
  private ci = new Float64Array(0);
  private dr = new Float64Array(0); // direct term Re(G)
  private f = new Float64Array(0);
  private lastGain: Float32Array<ArrayBufferLike> = new Float32Array(0);
  private exciters: Exciter[] = []; // due in this block, or already sounding
  private pending: Exciter[] = []; // waiting for a later frame
  private loadAcc = 0;
  private loadBlocks = 0;
  private lastLoadReport = 0;

  constructor(options?: { processorOptions?: { init?: InitMessage; hits?: HitMessage[] } }) {
    super();
    this.port.onmessage = (ev: MessageEvent<InMessage>) => this.onMessage(ev.data);
    // An offline render hands the whole performance over here, because a port message posted
    // before startRendering is not guaranteed to arrive before the first block.
    const opts = options?.processorOptions;
    if (opts?.init) this.init(opts.init);
    if (opts?.hits) for (const h of opts.hits) this.addHit(h);
  }

  private onMessage(msg: InMessage) {
    switch (msg.type) {
      case "init":
        this.init(msg);
        break;
      case "params":
        this.setParams(msg.freq, msg.sigma);
        break;
      case "gain":
        this.setGain(msg.gain);
        break;
      case "hit":
        this.addHit(msg);
        break;
      case "reset":
        this.zr.fill(0);
        this.zi.fill(0);
        this.exciters = [];
        this.pending = [];
        break;
    }
  }

  private init(msg: InitMessage) {
    const Nm = msg.freq.length;
    this.nModes = Nm;
    this.zr = new Float64Array(Nm);
    this.zi = new Float64Array(Nm);
    this.lamR = new Float64Array(Nm);
    this.lamI = new Float64Array(Nm);
    this.pr = new Float64Array(Nm);
    this.pi = new Float64Array(Nm);
    this.br = new Float64Array(Nm);
    this.bi = new Float64Array(Nm);
    this.cr = new Float64Array(Nm);
    this.ci = new Float64Array(Nm);
    this.dr = new Float64Array(Nm);
    this.f = new Float64Array(Nm);
    this.lastGain = msg.gain;
    this.setParams(msg.freq, msg.sigma);
    this.exciters = [];
    this.pending = [];
    this.warmUp();
  }

  /** Run realistic fake hits through the contact code path so that the JIT has compiled and
   * optimised it before the first real one; the first contact blocks were otherwise several
   * times over the deadline. */
  private warmUp() {
    const Nm = this.nModes;
    const scratch = new Float32Array(128);
    const fake = (P: number, kPoint: number, alpha: number, mass: number, velocity: number) => {
      const phi = new Float32Array(P * Nm), offsets = new Float32Array(P);
      for (let p = 0; p < P; p++) {
        for (let m = 0; m < Nm; m++) phi[p * Nm + m] = 0.7 * Math.sin(0.37 * m + 1.3 * p) * Math.cos(0.11 * m);
        offsets[p] = 1e-4 * (p - 0.5 * (P - 1));
      }
      this.addHit({ type: "hit", id: -1, phi, offsets, kPoint, alpha, mass, gap: 1e-4, velocity, silent: true });
    };
    fake(1, 1.8e8, 2.5, 0.3, 8);
    fake(7, 5e8 / 7, 1.5, 0.03, 7);
    for (let k = 0; k < 80; k++) {
      this.renderBlock(scratch, 0);
      if (k === 40) fake(7, 5e8 / 7, 1.5, 0.03, 5);
    }
    this.exciters = [];
    this.pending = [];
    this.zr.fill(0);
    this.zi.fill(0);
  }

  /** Poles and forcing coefficients of every mode; the ringing state is left alone. */
  private setParams(freq: Float64Array, sigma: Float64Array) {
    const Nm = this.nModes;
    if (freq.length !== Nm) return;
    const dt = 1 / sampleRate;
    for (let m = 0; m < Nm; m++) {
      const w0 = 2 * Math.PI * freq[m];
      const s = Math.min(sigma[m], 0.999 * w0);
      const wd = Math.sqrt(Math.max(w0 * w0 - s * s, 1e-6));
      const lr = -s, li = wd;
      this.lamR[m] = lr;
      this.lamI[m] = li;
      const e = Math.exp(lr * dt);
      const Pr = e * Math.cos(li * dt), Pi = e * Math.sin(li * dt);
      this.pr[m] = Pr;
      this.pi[m] = Pi;
      // b = (P - 1) / (lambda omega_d)
      const nr = Pr - 1, ni = Pi;
      const den = wd * (lr * lr + li * li);
      this.br[m] = (nr * lr + ni * li) / den;
      this.bi[m] = (ni * lr - nr * li) / den;
    }
    this.setGain(this.lastGain);
  }

  private setGain(gain: Float32Array) {
    const Nm = this.nModes;
    this.lastGain = gain;
    if (gain.length < 2 * Nm) return;
    for (let m = 0; m < Nm; m++) {
      const gr = gain[2 * m], gi = gain[2 * m + 1];
      const lr = this.lamR[m], li = this.lamI[m];
      const l2r = lr * lr - li * li, l2i = 2 * lr * li;
      this.cr[m] = gr * l2r - gi * l2i;
      this.ci[m] = gr * l2i + gi * l2r;
      this.dr[m] = gr;
    }
  }

  private addHit(msg: HitMessage) {
    const P = msg.offsets.length, Nm = this.nModes;
    if (Nm === 0 || msg.phi.length !== P * Nm) return;
    // one-step compliance c_pq = sum_m phi_m(p) phi_m(q) Im(b_m)
    const c = new Float64Array(P * P);
    const phi = msg.phi, bi = this.bi;
    for (let p = 0; p < P; p++) {
      for (let q = p; q < P; q++) {
        let s = 0;
        const op = p * Nm, oq = q * Nm;
        for (let m = 0; m < Nm; m++) s += bi[m] * phi[op + m] * phi[oq + m];
        c[p * P + q] = s;
        c[q * P + p] = s;
      }
    }
    const ex: Exciter = {
      id: msg.id, phi, offsets: msg.offsets, nPoints: P, kPoint: msg.kPoint, alpha: msg.alpha, mass: msg.mass,
      gap: msg.gap, velocity: msg.velocity, atFrame: msg.atFrame ?? 0, silent: msg.silent ?? false,
      released: false, zs: HOLD, vs: 0, startFrame: 0, peakForce: 0, contactSamples: 0, newtonIters: 0,
      touched: false, wp: new Float64Array(P), Fp: new Float64Array(P), c,
      J: new Float64Array(P * P), g: new Float64Array(P), dF: new Float64Array(P),
    };
    if (ex.atFrame > currentFrame + 128) this.pending.push(ex);
    else this.exciters.push(ex);
  }

  /** Surface height (with offsets) under the patch points from the current state; returns the max. */
  private surface(ex: Exciter): number {
    const Nm = this.nModes, zi = this.zi, phi = ex.phi, P = ex.nPoints, wp = ex.wp;
    let top = -Infinity;
    for (let p = 0; p < P; p++) {
      const o = p * Nm;
      let s = 0;
      for (let m = 0; m < Nm; m++) s += phi[o + m] * zi[m];
      s += ex.offsets[p];
      wp[p] = s;
      if (s > top) top = s;
    }
    return top;
  }

  /** Implicit contact forces for this sample given the free end-of-step surface heights ex.wp. */
  private solveContact(ex: Exciter, dt: number): number {
    const P = ex.nPoints, c = ex.c, F = ex.Fp, J = ex.J, g = ex.g, dF = ex.dF;
    const base = ex.zs + dt * ex.vs; // the exciter face at the end of the step, without contact
    const em = dt * dt / ex.mass;
    const k = ex.kPoint, a = ex.alpha;
    let touch = false;
    for (let p = 0; p < P; p++) if (ex.wp[p] > base) touch = true;
    if (!touch) {
      F.fill(0);
      return 0;
    }
    for (let it = 0; it < NEWTON_MAX; it++) {
      ex.newtonIters++;
      let Ftot = 0;
      for (let q = 0; q < P; q++) Ftot += F[q];
      let worst = 0;
      for (let p = 0; p < P; p++) {
        let d = ex.wp[p] - base - em * Ftot;
        for (let q = 0; q < P; q++) d -= c[p * P + q] * F[q];
        if (d > 0) {
          const dm1 = Math.pow(d, a - 1);
          g[p] = F[p] - k * dm1 * d;
          const s = k * a * dm1;
          for (let q = 0; q < P; q++) J[p * P + q] = (p === q ? 1 : 0) + s * (c[p * P + q] + em);
        } else {
          g[p] = F[p];
          for (let q = 0; q < P; q++) J[p * P + q] = p === q ? 1 : 0;
        }
      }
      solveSmall(J, g, dF, P);
      for (let p = 0; p < P; p++) {
        F[p] = Math.max(0, F[p] - dF[p]);
        worst = Math.max(worst, Math.abs(dF[p]));
      }
      let scale = 1;
      for (let p = 0; p < P; p++) scale = Math.max(scale, F[p]);
      if (worst < 1e-7 * scale) break;
    }
    let Ftot = 0;
    for (let p = 0; p < P; p++) Ftot += F[p];
    return Ftot;
  }

  process(_inputs: Float32Array[][], outputs: Float32Array[][]): boolean {
    const out = outputs[0][0];
    if (!out || this.nModes === 0) return true;
    const t0 = now();
    this.renderBlock(out, currentFrame);
    const cost = now() - t0;
    this.loadAcc += cost / (1000 * out.length / sampleRate);
    this.loadBlocks++;
    if (currentFrame - this.lastLoadReport >= sampleRate / 4) {
      this.port.postMessage({ type: "load", mean: this.loadAcc / this.loadBlocks, frame: currentFrame });
      this.loadAcc = 0;
      this.loadBlocks = 0;
      this.lastLoadReport = currentFrame;
    }
    return true;
  }

  private renderBlock(out: Float32Array, frame0: number): void {
    const B = out.length;
    const Nm = this.nModes;
    const zr = this.zr, zi = this.zi, cr = this.cr, ci = this.ci, pr = this.pr, pi = this.pi;
    const dt = 1 / sampleRate;

    // anything due in this block joins the working set; the rest costs nothing yet
    if (this.pending.length) {
      const end = frame0 + B;
      const still: Exciter[] = [];
      for (const ex of this.pending) {
        if (ex.atFrame < end) this.exciters.push(ex);
        else still.push(ex);
      }
      this.pending = still;
    }

    if (this.exciters.length === 0) {
      // free ringing: mode-outer over the block
      out.fill(0);
      for (let m = 0; m < Nm; m++) {
        let r = zr[m], i = zi[m];
        const a = pr[m], b = pi[m], c = cr[m], d = ci[m];
        for (let n = 0; n < B; n++) {
          const r2 = a * r - b * i;
          i = a * i + b * r;
          r = r2;
          out[n] += c * i + d * r;
        }
        zr[m] = r;
        zi[m] = i;
      }
    } else {
      const f = this.f, dr = this.dr, br = this.br, bi = this.bi;
      const exs = this.exciters;
      for (let n = 0; n < B; n++) {
        // free step of the bank
        for (let m = 0; m < Nm; m++) {
          const r = zr[m], i = zi[m];
          zr[m] = pr[m] * r - pi[m] * i;
          zi[m] = pr[m] * i + pi[m] * r;
        }
        // release whatever is due at this sample, just above the surface as it is now
        const frame = frame0 + n;
        for (let e = 0; e < exs.length; e++) {
          const ex = exs[e];
          if (ex.released || frame < ex.atFrame) continue;
          ex.zs = this.surface(ex) + ex.gap;
          ex.vs = -ex.velocity;
          ex.released = true;
          ex.startFrame = frame;
          if (!ex.silent) this.port.postMessage({ type: "hitStarted", id: ex.id, frame, time: frame / sampleRate });
        }
        // implicit contact forces, then the exciter update
        let any = false;
        for (let e = 0; e < exs.length; e++) {
          const ex = exs[e];
          if (!ex.released) continue;
          this.surface(ex);
          const F = this.solveContact(ex, dt);
          if (F > 0) {
            any = true;
            ex.touched = true;
            ex.contactSamples++;
            if (F > ex.peakForce) ex.peakForce = F;
          }
          ex.vs += dt * F / ex.mass;
          ex.zs += dt * ex.vs;
        }
        // modal forces, forced part of the step, output
        let y = 0;
        if (any) {
          f.fill(0);
          for (let e = 0; e < exs.length; e++) {
            const ex = exs[e], P = ex.nPoints, phi = ex.phi, F = ex.Fp;
            if (!ex.released) continue;
            for (let p = 0; p < P; p++) {
              const Fp = F[p];
              if (Fp === 0) continue;
              const o = p * Nm;
              for (let m = 0; m < Nm; m++) f[m] -= Fp * phi[o + m];
            }
          }
          for (let m = 0; m < Nm; m++) {
            const fm = f[m];
            const r = zr[m] + br[m] * fm, i = zi[m] + bi[m] * fm;
            zr[m] = r;
            zi[m] = i;
            y += cr[m] * i + ci[m] * r + dr[m] * fm;
          }
        } else {
          for (let m = 0; m < Nm; m++) y += cr[m] * zi[m] + ci[m] * zr[m];
        }
        out[n] = y;
      }
      // drop exciters that have rebounded clear of the surface
      const frameEnd = frame0 + B;
      this.exciters = exs.filter((ex) => {
        if (!ex.released) return true;
        const top = this.surface(ex);
        const done = ex.touched && ex.vs > 0 && ex.zs - top > RELEASE_CLEARANCE;
        const lost = !ex.touched && frameEnd - ex.startFrame > sampleRate * 0.1;
        if (done || lost) {
          if (!ex.silent) this.port.postMessage({
            type: "hitDone", id: ex.id, peakForce: ex.peakForce, contactTime: ex.contactSamples * dt,
            rebound: ex.vs, newtonIters: ex.newtonIters, contactSamples: ex.contactSamples,
          });
          return false;
        }
        return true;
      });
    }
  }
}

/** Solve A x = b for a small dense system (row-major A, n <= 16) by Gaussian elimination with
 * partial pivoting; A and b are destroyed. */
function solveSmall(A: Float64Array, b: Float64Array, x: Float64Array, n: number): void {
  for (let c = 0; c < n; c++) {
    let piv = c, best = Math.abs(A[c * n + c]);
    for (let r = c + 1; r < n; r++) {
      const v = Math.abs(A[r * n + c]);
      if (v > best) { best = v; piv = r; }
    }
    if (piv !== c) {
      for (let k = 0; k < n; k++) { const t = A[c * n + k]; A[c * n + k] = A[piv * n + k]; A[piv * n + k] = t; }
      const t = b[c]; b[c] = b[piv]; b[piv] = t;
    }
    const d = A[c * n + c] || 1e-300;
    for (let r = c + 1; r < n; r++) {
      const fct = A[r * n + c] / d;
      if (fct === 0) continue;
      for (let k = c; k < n; k++) A[r * n + k] -= fct * A[c * n + k];
      b[r] -= fct * b[c];
    }
  }
  for (let r = n - 1; r >= 0; r--) {
    let s = b[r];
    for (let k = r + 1; k < n; k++) s -= A[r * n + k] * x[k];
    x[r] = s / (A[r * n + r] || 1e-300);
  }
}

registerProcessor("cymbal-modal", CymbalModalProcessor);