// MQ Synthesizer
// Replica oscillator bank for sinusoidal synthesis, plus two-pole resonator mode

// Evaluate cubic bezier curve at time t
function evalBezier(curve, t) {
  const dt = curve.t3 - curve.t0;
  if (dt <= 0) return curve.v0;
  let u = (t - curve.t0) / dt;
  u = Math.max(0, Math.min(1, u));
  const u1 = 1.0 - u;
  return u1*u1*u1 * curve.v0 +
         3*u1*u1*u * curve.v1 +
         3*u1*u*u  * curve.v2 +
         u*u*u     * curve.v3;
}

// Deterministic LCG PRNG
function randFloat(seed, min, max) {
  seed = (1664525 * seed + 1013904223) % 0x100000000;
  return min + (seed / 0x100000000) * (max - min);
}

// Synthesize audio from MQ partials
// partials: array of {freqCurve, ampCurve, replicas?, resonator?}
// replicas:  {offsets, decay_alpha, jitter, spread_above, spread_below}
// resonator: {enabled, r, gainComp}  — two-pole resonator mode per partial
// integratePhase: true = accumulate 2π*f/SR per sample (correct for varying freq)
//                 false = 2π*f*t (simpler, only correct for constant freq)
// options.k1: LP coefficient in (0,1] — omit to bypass
// options.k2: HP coefficient in (0,1] — omit to bypass
function synthesizeMQ(partials, sampleRate, duration, integratePhase = true, options = {}) {
  const numSamples = Math.floor(sampleRate * duration);
  const pcm = new Float32Array(numSamples);

  const jitterMult = options.disableJitter ? 0 : 1;
  const spreadMult = options.disableSpread ? 0 : 1;

  const defaultReplicas = {
    offsets:       [1.0],
    decay_alpha:   0.1,
    jitter:        0.05,
    spread_above:  0.02,
    spread_below:  0.02
  };

  // Pre-build per-partial configs with fixed spread/jitter and phase accumulators
  const configs = [];
  for (let p = 0; p < partials.length; ++p) {
    const partial = partials[p];
    const fc = partial.freqCurve;
    const ac = partial.ampCurve;

    if ((partial.resonator && partial.resonator.enabled) || options.forceResonator) {
      // --- Two-pole resonator mode ---
      // Driven by band-limited noise scaled by amp curve.
      // r controls pole radius (bandwidth): r→1 = narrow, r→0 = wide.
      // gainNorm = sqrt(1 - r²) normalises steady-state output power to ~A.
      const res      = partial.resonator || {};
      const r        = options.forceRGain ? Math.min(0.9999, Math.max(0, options.globalR))
                     : (res.r        != null ? Math.min(0.9999, Math.max(0, res.r)) : 0.995);
      const gainComp = options.forceRGain ? options.globalGain
                     : (res.gainComp != null ? res.gainComp : 1.0);
      const gainNorm = Math.sqrt(Math.max(0, 1.0 - r * r));
      configs.push({
        mode: 'resonator',
        fc, ac,
        r, gainComp, gainNorm,
        y1: 0.0, y2: 0.0,
        noiseSeed: ((p * 1664525 + 1013904223) & 0xFFFFFFFF) >>> 0
      });
    } else {
      // --- Sinusoidal (replica) mode ---
      const rep = partial.replicas != null ? partial.replicas : defaultReplicas;
      const offsets      = rep.offsets      != null ? rep.offsets      : [1.0];
      const decay_alpha  = rep.decay_alpha  != null ? rep.decay_alpha  : 0.0;
      const jitter       = rep.jitter       != null ? rep.jitter       : 0.0;
      const spread_above = rep.spread_above != null ? rep.spread_above : 0.0;
      const spread_below = rep.spread_below != null ? rep.spread_below : 0.0;

      const replicaData = [];
      for (let r = 0; r < offsets.length; ++r) {
        const spread    = spreadMult * randFloat(p * 67890 + r * 999, -spread_below, spread_above);
        const initPhase = randFloat(p * 67890 + r, 0.0, 1.0) * (jitter * jitterMult) * 2.0 * Math.PI;
        replicaData.push({ratio: offsets[r], spread, phase: initPhase});
      }
      configs.push({ mode: 'sinusoid', fc, ac, decay_alpha, replicaData });
    }
  }

  for (let i = 0; i < numSamples; ++i) {
    const t = i / sampleRate;
    let sample = 0.0;

    for (let p = 0; p < configs.length; ++p) {
      const cfg = configs[p];
      const {fc, ac} = cfg;

      if (cfg.mode === 'resonator') {
        if (t < fc.t0 || t > fc.t3) { cfg.y1 = 0.0; cfg.y2 = 0.0; continue; }

        const f0    = evalBezier(fc, t);
        const A     = evalBezier(ac, t);
        const omega = 2.0 * Math.PI * f0 / sampleRate;
        const b1    = 2.0 * cfg.r * Math.cos(omega);

        // LCG noise excitation (deterministic per-partial)
        cfg.noiseSeed = (Math.imul(1664525, cfg.noiseSeed) + 1013904223) >>> 0;
        const noise   = cfg.noiseSeed / 0x100000000 * 2.0 - 1.0;

        const x  = A * cfg.gainNorm * noise;
        const y  = b1 * cfg.y1 - cfg.r * cfg.r * cfg.y2 + x;
        cfg.y2   = cfg.y1;
        cfg.y1   = y;
        sample  += y * cfg.gainComp;

      } else {
        if (t < fc.t0 || t > fc.t3) continue;

        const f0 = evalBezier(fc, t);
        const A0 = evalBezier(ac, t);
        const {decay_alpha, replicaData} = cfg;

        for (let r = 0; r < replicaData.length; ++r) {
          const rep = replicaData[r];
          const f = f0 * rep.ratio * (1.0 + rep.spread);
          const A = A0 * Math.exp(-decay_alpha * Math.abs(f - f0));

          let phase;
          if (integratePhase) {
            rep.phase += 2.0 * Math.PI * f / sampleRate;
            phase = rep.phase;
          } else {
            phase = 2.0 * Math.PI * f * t + rep.phase;
          }

          sample += A * Math.sin(phase);
        }
      }
    }

    pcm[i] = sample;
  }

  // Post-synthesis filters (applied before normalization)
  // LP: y[n] = k1*x[n] + (1-k1)*y[n-1]  — options.k1 in (0,1], omit to bypass
  // HP: y[n] = k2*(y[n-1] + x[n] - x[n-1]) — options.k2 in (0,1], omit to bypass
  if (options.k1 != null) {
    const k1 = Math.max(0, Math.min(1, options.k1));
    let y = 0.0;
    for (let i = 0; i < numSamples; ++i) {
      y = k1 * pcm[i] + (1.0 - k1) * y;
      pcm[i] = y;
    }
  }
  if (options.k2 != null) {
    const k2 = Math.max(0, Math.min(1, options.k2));
    let y = 0.0, xp = 0.0;
    for (let i = 0; i < numSamples; ++i) {
      const x = pcm[i];
      y = k2 * (y + x - xp);
      xp = x;
      pcm[i] = y;
    }
  }

  // Normalize
  let maxAbs = 0;
  for (let i = 0; i < numSamples; ++i) maxAbs = Math.max(maxAbs, Math.abs(pcm[i]));
  if (maxAbs > 1.0) {
    for (let i = 0; i < numSamples; ++i) pcm[i] /= maxAbs;
  }

  return pcm;
}