import type { Anchor, BeadSystem, Excitation } from '../../sim/types';
import { tridiagEigenvalues } from '../../sim/tridiag';

export type EndType = 'clamped' | 'sliding' | 'free';
export type ExcitationType = 'pluck' | 'strike' | 'hammer';

export interface StringParams {
  // geometry
  n: number; // beads
  spacing: number; // m, bead spacing at equilibrium (also bead-to-wall)
  // properties
  mass: number; // kg per bead
  k: number; // spring stiffness, N/m
  tension: number; // N, equilibrium tension (sets rest length L0 = spacing - T/k)
  leftEnd: EndType;
  rightEnd: EndType;
  pins: string; // 1-based bead numbers, comma separated
  // damping (0 = off)
  t60Low: number; // s
  fLow: number; // Hz
  t60High: number; // s
  fHigh: number; // Hz
  // excitation
  excitation: ExcitationType;
  excBead: number; // 1-based
  excWidth: number; // beads (Gaussian sigma); 0 = one bead
  pluckAmp: number; // m
  strikeVel: number; // m/s
  hammerForce: number; // N, peak total
  hammerMs: number; // ms
  excAngle: number; // degrees from y toward z
  // pickup
  pickupMode: 'bead' | 'sum';
  pickupBead: number; // 1-based
  pickupAxis: 'y' | 'z' | 'x';
  // simulation
  duration: number;
  fs: number;
  cfl: number;
  frameRate: number;
}

export function parsePins(s: string, n: number): number[] {
  return s
    .split(/[\s,;]+/)
    .filter((t) => t.length > 0)
    .map((t) => Math.round(Number(t)) - 1)
    .filter((i) => Number.isFinite(i) && i >= 0 && i < n);
}

export interface Derived {
  restLength: number;
  prestretch: number; // (L - L0) / L0
  length: number; // wall to wall (m)
  linearDensity: number; // kg/m
  waveSpeedT: number; // transverse, continuum limit (m/s)
  waveSpeedL: number; // longitudinal, continuum limit (m/s)
  transverseHz: number[]; // linear modes (lowest first)
  longitudinalHz: number[];
  sigma0: number;
  sigma1: number;
  dashpot: number;
  warnings: string[];
}

// Linear damping sigma(omega) = sigma0 + sigma1 omega^2 for transverse modes,
// fitted to two decay times T60 = 3 ln(10) / sigma.
export function dampingFromT60(
  p: Pick<StringParams, 't60Low' | 'fLow' | 't60High' | 'fHigh'>,
): { sigma0: number; sigma1: number; note?: string } {
  if (!(p.t60Low > 0) && !(p.t60High > 0)) return { sigma0: 0, sigma1: 0 };
  const ln1000 = 3 * Math.LN10;
  if (!(p.t60High > 0) || !(p.fHigh > p.fLow)) return { sigma0: ln1000 / p.t60Low, sigma1: 0 };
  if (!(p.t60Low > 0)) {
    const w = 2 * Math.PI * p.fHigh;
    return { sigma0: 0, sigma1: ln1000 / p.t60High / (w * w) };
  }
  const s1 = ln1000 / p.t60Low;
  const s2 = ln1000 / p.t60High;
  const w1 = 2 * Math.PI * p.fLow;
  const w2 = 2 * Math.PI * p.fHigh;
  let sigma1 = (s2 - s1) / (w2 * w2 - w1 * w1);
  let sigma0 = s1 - sigma1 * w1 * w1;
  let note: string | undefined;
  if (sigma1 < 0) {
    sigma1 = 0;
    sigma0 = s1;
    note = 'High-frequency T60 is longer than low-frequency T60; using frequency-independent damping.';
  } else if (sigma0 < 0) {
    sigma0 = 0;
    sigma1 = s1 / (w1 * w1);
    note = 'Decay times imply negative mass damping; clamped to zero.';
  }
  return { sigma0, sigma1, note };
}

// Linearized modes about the straight equilibrium. Transverse: each spring acts
// with stiffness T/L (zero to a sliding ring, which follows the bead).
// Longitudinal: stiffness k (to clamped or sliding walls; nothing to a free end).
function linearModes(p: StringParams, count: number): { trans: number[]; long: number[] } {
  const pinned = new Set(parsePins(p.pins, p.n));
  const idx: number[] = [];
  for (let i = 0; i < p.n; i++) if (!pinned.has(i)) idx.push(i);
  const build = (kSpring: number, wallLeft: number, wallRight: number) => {
    const d = new Float64Array(idx.length);
    const e = new Float64Array(Math.max(0, idx.length - 1));
    idx.forEach((i, r) => {
      let diag = 0;
      diag += i === 0 ? wallLeft : kSpring;
      diag += i === p.n - 1 ? wallRight : kSpring;
      d[r] = diag / p.mass;
      if (r < idx.length - 1 && idx[r + 1] === i + 1) e[r] = -kSpring / p.mass;
    });
    return tridiagEigenvalues(d, e, count)
      .filter((l) => l > 1e-9 * Math.max(...d, 1))
      .map((l) => Math.sqrt(l) / (2 * Math.PI));
  };
  const kt = p.tension / p.spacing;
  const tWall = (end: EndType) => (end === 'clamped' ? kt : 0);
  const lWall = (end: EndType) => (end === 'free' ? 0 : p.k);
  return {
    trans: build(kt, tWall(p.leftEnd), tWall(p.rightEnd)),
    long: build(p.k, lWall(p.leftEnd), lWall(p.rightEnd)),
  };
}

export function derive(p: StringParams): Derived {
  const warnings: string[] = [];
  const restLength = p.spacing - p.tension / p.k;
  if (!(restLength > 0)) warnings.push('Tension too high for this stiffness: rest length L0 = L - T/k is not positive.');
  if (p.tension > 0 && (p.leftEnd === 'free' || p.rightEnd === 'free'))
    warnings.push('A free end cannot hold tension: the string is not in equilibrium and will contract. Use sliding ends (tension kept, free transverse motion) or set T = 0.');
  if (p.tension < 0) warnings.push('Negative tension (compression): the straight string is unstable to buckling.');
  const linearDensity = p.mass / p.spacing;
  const { sigma0, sigma1, note } = dampingFromT60(p);
  if (note) warnings.push(note);
  const kt = p.tension > 0 ? p.tension / p.spacing : p.k;
  const modes = linearModes(p, 24);
  return {
    restLength,
    prestretch: p.spacing / restLength - 1,
    length: (p.n + 1) * p.spacing,
    linearDensity,
    waveSpeedT: Math.sqrt(Math.max(p.tension, 0) / linearDensity),
    waveSpeedL: Math.sqrt((p.k * p.spacing) / linearDensity),
    transverseHz: modes.trans,
    longitudinalHz: modes.long,
    sigma0,
    sigma1,
    dashpot: 2 * sigma1 * kt,
    warnings,
  };
}

export function buildSystem(p: StringParams, d: Derived): BeadSystem {
  const n = p.n;
  const rest = new Float64Array(3 * n);
  for (let i = 0; i < n; i++) rest[3 * i] = (i + 1) * p.spacing;
  const mass = new Float64Array(n).fill(p.mass);
  const pinned = new Uint8Array(n);
  for (const i of parsePins(p.pins, n)) pinned[i] = 1;
  const anchors: Anchor[] = [];
  const si: number[] = [];
  const sj: number[] = [];
  const addEnd = (end: EndType, x: number, bead: number) => {
    if (end === 'free') return;
    anchors.push({ pos: [x, 0, 0], slideBead: end === 'sliding' ? bead : -1 });
    const node = n + anchors.length - 1;
    // anchor on the left of bead 0 is si; keep orientation left -> right
    if (bead === 0) {
      si.push(node);
      sj.push(0);
    } else {
      si.push(bead);
      sj.push(node);
    }
  };
  addEnd(p.leftEnd, 0, 0);
  for (let i = 0; i < n - 1; i++) {
    si.push(i);
    sj.push(i + 1);
  }
  addEnd(p.rightEnd, (n + 1) * p.spacing, n - 1);
  const ns = si.length;
  return {
    n,
    rest,
    mass,
    pinned,
    anchors,
    si: Int32Array.from(si),
    sj: Int32Array.from(sj),
    sk: new Float64Array(ns).fill(p.k),
    sL0: new Float64Array(ns).fill(d.restLength),
    sc: new Float64Array(ns).fill(d.dashpot),
    sigma0: d.sigma0,
  };
}

function gaussianWeights(n: number, center: number, width: number): Float64Array {
  const w = new Float64Array(n);
  if (!(width > 0)) {
    w[center] = 1;
    return w;
  }
  for (let i = 0; i < n; i++) w[i] = Math.exp(-0.5 * ((i - center) / width) ** 2);
  return w;
}

const clampIndex = (i: number, n: number) => Math.min(n - 1, Math.max(0, Math.round(i) - 1));

export function buildExcitation(p: StringParams): Excitation {
  const n = p.n;
  const u0 = new Float64Array(3 * n);
  const v0 = new Float64Array(3 * n);
  const th = (p.excAngle * Math.PI) / 180;
  const dy = Math.cos(th);
  const dz = Math.sin(th);
  const c = clampIndex(p.excBead, n);
  if (p.excitation === 'pluck') {
    // triangle from the nearest support on each side (pinned bead or wall) to the plucked bead
    const pins = new Set(parsePins(p.pins, n));
    let left = -1;
    for (let i = c - 1; i >= 0; i--) if (pins.has(i)) { left = i; break; }
    let right = n;
    for (let i = c + 1; i < n; i++) if (pins.has(i)) { right = i; break; }
    for (let i = left + 1; i < right; i++) {
      const s = i <= c ? (i - left) / (c - left) : (right - i) / (right - c);
      u0[3 * i + 1] = p.pluckAmp * s * dy;
      u0[3 * i + 2] = p.pluckAmp * s * dz;
    }
    return { u0, v0 };
  }
  const w = gaussianWeights(n, c, p.excWidth);
  if (p.excitation === 'strike') {
    for (let i = 0; i < n; i++) {
      v0[3 * i + 1] = p.strikeVel * w[i] * dy;
      v0[3 * i + 2] = p.strikeVel * w[i] * dz;
    }
    return { u0, v0 };
  }
  let sum = 0;
  for (let i = 0; i < n; i++) sum += w[i];
  const shape = new Float64Array(3 * n);
  for (let i = 0; i < n; i++) {
    shape[3 * i + 1] = (w[i] / sum) * dy;
    shape[3 * i + 2] = (w[i] / sum) * dz;
  }
  return { u0, v0, force: { shape, F0: p.hammerForce, tau: p.hammerMs / 1000 } };
}

export function buildPickup(p: StringParams): Float64Array {
  const w = new Float64Array(3 * p.n);
  const axis = p.pickupAxis === 'x' ? 0 : p.pickupAxis === 'y' ? 1 : 2;
  if (p.pickupMode === 'sum') for (let i = 0; i < p.n; i++) w[3 * i + axis] = 1;
  else w[3 * clampIndex(p.pickupBead, p.n) + axis] = 1;
  return w;
}