import type { BeadSystem, Bending, Excitation } from '../../sim/types';
import { hingeAngleGrad } from '../../sim/hinge';
import { dampingFromT60 } from '../bead-string/model';

// A cymbal: a thin bronze plate, round, free at the rim, with a hole in the
// centre. Nothing holds it: it is simulated as a free plate whose rigid-body
// motion is removed (see BeadSystem.freeRigid), so it rings with no support
// stiffness or damping. Beads on a triangular lattice in
// the x-z plane, lifted to the profile height y(r) (flat by default; optional
// bow and bell). In-plane springs on every edge; plate bending from the
// dihedral angles across interior edges, per triangle (see Bending in sim/types). Vibration is mostly along y.

export type CymbalExcitation = 'stick' | 'mallet' | 'strike' | 'force';

export interface CymbalParams {
  // geometry
  diameter: number; // m
  rings: number; // beads from the centre to the rim
  holeRadius: number; // m
  bowHeight: number; // m, dome rise from rim to centre (0 = flat)
  bellDiameter: number; // m
  bellHeight: number; // m (0 = no bell)
  // material
  thickness: number; // m
  youngs: number; // Pa
  density: number; // kg/m^3
  soften: number; // in-plane stiffness divided by this (1 = real bronze)
  // damping
  t60Low: number;
  fLow: number;
  t60High: number;
  fHigh: number;
  // excitation
  excitation: CymbalExcitation;
  excR: number; // fraction of radius
  excAngle: number; // degrees, azimuth in the x-z plane
  hitAngle: number; // degrees from the surface normal, tilted toward the rim
  excWidth: number; // m, Gaussian sigma of the patch (0 = one bead)
  stickMass: number;
  stickSpeed: number;
  stickK: number; // N/m^1.5 (Hertz)
  malletMass: number;
  malletSpeed: number;
  malletK: number; // N/m^p
  malletP: number;
  friction: number;
  strikeVel: number;
  forcePeak: number;
  forceMs: number;
  // microphone
  micDistance: number; // m
  micElevation: number; // degrees above the plane
  micAzimuth: number; // degrees
  // simulation
  duration: number;
  fs: number;
  cfl: number;
  frameRate: number;
}

export function profileHeight(p: CymbalParams, r: number): number {
  const R = p.diameter / 2;
  let y = p.bowHeight * (1 - (r / R) ** 2);
  const rb = p.bellDiameter / 2;
  if (p.bellHeight > 0 && rb > 0 && r < rb) y += p.bellHeight * 0.5 * (1 + Math.cos((Math.PI * r) / rb));
  return y;
}

export interface CymbalMesh {
  a: number;
  pos: [number, number][]; // (x, z)
  tris: [number, number, number][];
  edges: { i: number; j: number; opp: number[] }[]; // opposite vertices of the 1 or 2 triangles
}

export function cymbalMesh(R: number, rings: number, hole: number): CymbalMesh {
  const a = R / (rings + 0.5);
  const at = (i: number, j: number): [number, number] => [a * (i + 0.5 * j), a * (Math.sqrt(3) / 2) * j];
  const index = new Map<string, number>();
  const pos: [number, number][] = [];
  const M = rings + 2;
  for (let j = -M; j <= M; j++)
    for (let i = -M - Math.abs(j); i <= M + Math.abs(j); i++) {
      const r = Math.hypot(...at(i, j));
      if (r <= R - 0.5 * a + 1e-9 && !(r < hole)) {
        index.set(`${i},${j}`, pos.length);
        pos.push(at(i, j));
      }
    }
  const id = (i: number, j: number) => index.get(`${i},${j}`);
  const tris: [number, number, number][] = [];
  for (const k of index.keys()) {
    const [i, j] = k.split(',').map(Number);
    const t1 = [id(i, j), id(i + 1, j), id(i, j + 1)];
    const t2 = [id(i + 1, j), id(i + 1, j + 1), id(i, j + 1)];
    for (const t of [t1, t2]) if (t.every((v) => v !== undefined)) tris.push(t as [number, number, number]);
  }
  const em = new Map<string, { i: number; j: number; opp: number[] }>();
  for (const [x, y, z] of tris)
    for (const [u, v, o] of [[x, y, z], [y, z, x], [z, x, y]]) {
      const key = u < v ? `${u},${v}` : `${v},${u}`;
      let e = em.get(key);
      if (!e) em.set(key, (e = { i: Math.min(u, v), j: Math.max(u, v), opp: [] }));
      e.opp.push(o);
    }
  // scale so the triangulated area equals the disk's (minus the hole): the
  // outermost row of beads otherwise leaves the plate about a/2 short
  let meshArea = 0;
  for (const [u, v, w] of tris) {
    const [p, q, r] = [pos[u], pos[v], pos[w]];
    meshArea += 0.5 * Math.abs((q[0] - p[0]) * (r[1] - p[1]) - (r[0] - p[0]) * (q[1] - p[1]));
  }
  const f = Math.sqrt((Math.PI * (R * R - hole * hole)) / meshArea);
  for (const p of pos) {
    p[0] *= f;
    p[1] *= f;
  }
  return { a: a * f, pos, tris, edges: [...em.values()] };
}

export interface CymbalDerived {
  mesh: CymbalMesh;
  n: number;
  rest: Float64Array; // 3n
  area: Float64Array; // n, lumped area per bead
  normal: Float64Array; // 3n, upward unit normals
  mass: Float64Array;
  k: number; // interior spring stiffness
  D: number; // bending stiffness
  sigma0: number;
  sigma1: number;
  waveSpeedP: number; // in-plane P waves (after softening)
  bendingCoef: number; // sqrt(D / rho t): bending waves omega = coef * k^2
  warnings: string[];
}

// Triangular lattice, spacing a, Poisson ratio 1/3:
//   in-plane: spring stiffness k = sqrt(3)/2 E t (an edge with one triangle gets half)
//   bending:  D = E t^3 / (12 (1 - nu^2)), through per-triangle shape operators
export function deriveCymbal(p: CymbalParams): CymbalDerived {
  const R = p.diameter / 2;
  const mesh = cymbalMesh(R, Math.max(2, Math.round(p.rings)), p.holeRadius);
  const n = mesh.pos.length;
  const rest = new Float64Array(3 * n);
  mesh.pos.forEach(([x, z], i) => {
    rest[3 * i] = x;
    rest[3 * i + 1] = profileHeight(p, Math.hypot(x, z));
    rest[3 * i + 2] = z;
  });
  const area = new Float64Array(n);
  const normal = new Float64Array(3 * n);
  for (const [u, v, w] of mesh.tris) {
    const e1 = [0, 1, 2].map((c) => rest[3 * v + c] - rest[3 * u + c]);
    const e2 = [0, 1, 2].map((c) => rest[3 * w + c] - rest[3 * u + c]);
    let nx = e1[1] * e2[2] - e1[2] * e2[1];
    let ny = e1[2] * e2[0] - e1[0] * e2[2];
    let nz = e1[0] * e2[1] - e1[1] * e2[0];
    const A2 = Math.hypot(nx, ny, nz);
    if (ny < 0) [nx, ny, nz] = [-nx, -ny, -nz];
    for (const b of [u, v, w]) {
      area[b] += A2 / 6;
      normal[3 * b] += nx;
      normal[3 * b + 1] += ny;
      normal[3 * b + 2] += nz;
    }
  }
  for (let i = 0; i < n; i++) {
    const l = Math.hypot(normal[3 * i], normal[3 * i + 1], normal[3 * i + 2]) || 1;
    for (let c = 0; c < 3; c++) normal[3 * i + c] /= l;
  }
  const rhoT = p.density * p.thickness;
  const mass = area.map((A) => rhoT * A);
  const Et = (p.youngs * p.thickness) / Math.max(p.soften, 1e-9);
  const k = (Math.sqrt(3) / 2) * Et;
  const nu = 1 / 3;
  const D = (p.youngs * p.thickness ** 3) / (12 * (1 - nu * nu));
  const warnings: string[] = [];
  if (p.soften > 1)
    warnings.push(
      `In-plane stiffness is softened ${p.soften}×. A flat plate bends the same at small amplitude; curved shapes lose curvature stiffening and big hits exaggerate the nonlinear coupling.`,
    );
  const { sigma0, sigma1, note } = dampingFromT60(p);
  if (note) warnings.push(note);
  const lambdaPlus2Mu = (3 * Math.sqrt(3) * k) / 4;
  return {
    mesh,
    n,
    rest,
    area,
    normal,
    mass,
    k,
    D,
    sigma0,
    sigma1,
    waveSpeedP: Math.sqrt(lambdaPlus2Mu / rhoT),
    bendingCoef: Math.sqrt(D / rhoT),
    warnings,
  };
}

// highest frequency whose bending wavelength spans six lattice spacings (where
// the lattice is still a fair plate)
export function bendingLimitHz(d: CymbalDerived): number {
  const kb = (2 * Math.PI) / (6 * d.mesh.a);
  return (d.bendingCoef * kb * kb) / (2 * Math.PI);
}

export function buildCymbal(d: CymbalDerived): BeadSystem {
  const { n, rest, mesh } = d;
  const si: number[] = [];
  const sj: number[] = [];
  const sk: number[] = [];
  const sL0: number[] = [];
  const sc: number[] = [];
  const len = (i: number, j: number) =>
    Math.hypot(rest[3 * j] - rest[3 * i], rest[3 * j + 1] - rest[3 * i + 1], rest[3 * j + 2] - rest[3 * i + 2]);
  for (const e of mesh.edges) {
    si.push(e.i);
    sj.push(e.j);
    sk.push((d.k / 2) * e.opp.length);
    sL0.push(len(e.i, e.j));
    sc.push(0);
  }
  const bending = buildBending(d, 1 / 3);
  return {
    n,
    rest,
    mass: d.mass,
    pinned: new Uint8Array(n),
    anchors: [],
    freeRigid: true,
    si: Int32Array.from(si),
    sj: Int32Array.from(sj),
    sk: Float64Array.from(sk),
    sL0: Float64Array.from(sL0),
    sc: Float64Array.from(sc),
    sigma0: d.sigma0,
    bending,
  };
}

export function hitPoint(p: CymbalParams): [number, number] {
  const r = (p.diameter / 2) * p.excR;
  const th = (p.excAngle * Math.PI) / 180;
  return [r * Math.cos(th), r * Math.sin(th)];
}

function nearest(d: CymbalDerived, x: number, z: number): number {
  let best = 0;
  let bd = Infinity;
  d.mesh.pos.forEach(([px, pz], i) => {
    const q = (px - x) ** 2 + (pz - z) ** 2;
    if (q < bd) [bd, best] = [q, i];
  });
  return best;
}

// patch weights (sum 1), the inward normal there and the hit direction
// (tilted from the normal toward the rim by hitAngle)
export function hitGeometry(p: CymbalParams, d: CymbalDerived) {
  const [cx, cz] = hitPoint(p);
  const w = new Float64Array(d.n);
  if (p.excWidth > 0)
    d.mesh.pos.forEach(([x, z], i) => {
      const r2 = ((x - cx) ** 2 + (z - cz) ** 2) / (p.excWidth * p.excWidth);
      w[i] = r2 < 9 ? Math.exp(-0.5 * r2) : 0;
    });
  if (!w.some((v) => v > 0)) w[nearest(d, cx, cz)] = 1;
  const sum = w.reduce((s, v) => s + v, 0);
  for (let i = 0; i < d.n; i++) w[i] /= sum;
  const up = [0, 0, 0];
  for (let i = 0; i < d.n; i++) for (let c = 0; c < 3; c++) up[c] += w[i] * d.normal[3 * i + c];
  const ul = Math.hypot(...up);
  const N = up.map((v) => v / ul);
  // outward tangent: radial direction minus its normal part
  const rl = Math.hypot(cx, cz) || 1;
  const er = [cx / rl, 0, cz / rl];
  const dot = er[0] * N[0] + er[2] * N[2];
  let t = [er[0] - dot * N[0], -dot * N[1], er[2] - dot * N[2]];
  const tl = Math.hypot(...t) || 1;
  t = t.map((v) => v / tl);
  const al = (p.hitAngle * Math.PI) / 180;
  const inward: [number, number, number] = [-N[0], -N[1], -N[2]];
  const dir: [number, number, number] = [0, 1, 2].map((c) => Math.cos(al) * inward[c] + Math.sin(al) * t[c]) as [number, number, number];
  return { w, inward, dir };
}

export function buildCymbalExcitation(p: CymbalParams, d: CymbalDerived): Excitation {
  const n = d.n;
  const u0 = new Float64Array(3 * n);
  const v0 = new Float64Array(3 * n);
  const { w, inward, dir } = hitGeometry(p, d);
  const wmax = Math.max(...w);
  if (p.excitation === 'strike') {
    for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) v0[3 * i + c] = (p.strikeVel * w[i] * dir[c]) / wmax;
    return { u0, v0 };
  }
  if (p.excitation === 'force') {
    const shape = new Float64Array(3 * n);
    for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) shape[3 * i + c] = w[i] * dir[c];
    return { u0, v0, force: { shape, F0: p.forcePeak, tau: p.forceMs / 1000 } };
  }
  const stick = p.excitation === 'stick';
  const speed = stick ? p.stickSpeed : p.malletSpeed;
  return {
    u0,
    v0,
    mallet: {
      weights: w,
      normal: inward,
      vel0: [dir[0] * speed, dir[1] * speed, dir[2] * speed],
      mass: stick ? p.stickMass : p.malletMass,
      K: stick ? p.stickK : p.malletK,
      p: stick ? 1.5 : p.malletP,
      mu: p.friction,
      reach: Math.max(2.5 * p.excWidth, 0.75 * d.mesh.a),
    },
  };
}

const RHO_AIR = 1.2;
const C_AIR = 343;

export function micDirection(p: CymbalParams): [number, number, number] {
  const el = (p.micElevation * Math.PI) / 180;
  const az = (p.micAzimuth * Math.PI) / 180;
  return [Math.cos(el) * Math.cos(az), Math.sin(el), Math.cos(el) * Math.sin(az)];
}

// Rayleigh integral: p(t) = rho / (2 pi) sum_i A_i a_n,i(t - R_i / c) / R_i,
// with a_n the acceleration along the normal on the microphone's side.
export function buildMic(p: CymbalParams, d: CymbalDerived): { weights: Float64Array; delay: Float64Array } {
  const [mx, my, mz] = micDirection(p).map((v) => v * p.micDistance);
  const side = p.micElevation >= 0 ? 1 : -1;
  const weights = new Float64Array(3 * d.n);
  const delay = new Float64Array(d.n);
  for (let i = 0; i < d.n; i++) {
    const R = Math.hypot(mx - d.rest[3 * i], my - d.rest[3 * i + 1], mz - d.rest[3 * i + 2]);
    delay[i] = (R / C_AIR) * p.fs;
    const g = (side * RHO_AIR * d.area[i]) / (2 * Math.PI * R);
    for (let c = 0; c < 3; c++) weights[3 * i + c] = g * d.normal[3 * i + c];
  }
  return { weights, delay };
}

// Hinges on interior edges, oriented so a positive change folds the surface
// toward its upward normal; per-triangle plate-energy matrices, condensed over
// free edges (the missing angle takes the value that leaves no moment there).
function buildBending(d: CymbalDerived, nu: number): Bending {
  const { rest, mesh, normal } = d;
  const key = (i: number, j: number) => (i < j ? `${i},${j}` : `${j},${i}`);
  const hingeOf = new Map<string, number>();
  const w0: number[] = [], w1: number[] = [], e0: number[] = [], e1: number[] = [];
  for (const e of mesh.edges)
    if (e.opp.length === 2) {
      hingeOf.set(key(e.i, e.j), w0.length);
      w0.push(e.opp[0]);
      w1.push(e.opp[1]);
      e0.push(e.i);
      e1.push(e.j);
    }
  const nh = w0.length;
  const g = new Float64Array(13);
  const theta0 = new Float64Array(nh);
  const sign = new Float64Array(nh);
  const lift = 1e-4 * mesh.a;
  for (let e = 0; e < nh; e++) {
    theta0[e] = hingeAngleGrad(rest, w0[e], w1[e], e0[e], e1[e], g);
    const Y = rest.slice();
    for (const w of [w0[e], w1[e]]) for (let c = 0; c < 3; c++) Y[3 * w + c] += lift * normal[3 * w + c];
    sign[e] = Math.sign(hingeAngleGrad(Y, w0[e], w1[e], e0[e], e1[e], g) - theta0[e]) || 1;
  }
  const nt = mesh.tris.length;
  const tri = new Int32Array(3 * nt).fill(-1);
  const M = new Float64Array(9 * nt);
  const P = (i: number) => [rest[3 * i], rest[3 * i + 1], rest[3 * i + 2]];
  mesh.tris.forEach(([u, v, w], t) => {
    const ed = [[v, w], [w, u], [u, v]];
    const vec = ed.map(([i, j]) => P(j).map((x, c) => x - P(i)[c]));
    const len = vec.map((q) => Math.hypot(...q));
    const tu = vec.map((q, k) => q.map((x) => x / len[k]));
    const [x1, x2] = [vec[2], vec[1].map((x) => -x)]; // u->v and u->w
    const A = 0.5 * Math.hypot(x1[1] * x2[2] - x1[2] * x2[1], x1[2] * x2[0] - x1[0] * x2[2], x1[0] * x2[1] - x1[1] * x2[0]);
    const beta = len.map((l) => l / (2 * A));
    const full = new Float64Array(9);
    for (let a = 0; a < 3; a++)
      for (let b = 0; b < 3; b++) {
        const c2 = (tu[a][0] * tu[b][0] + tu[a][1] * tu[b][1] + tu[a][2] * tu[b][2]) ** 2;
        full[3 * a + b] = A * d.D * beta[a] * beta[b] * (nu + (1 - nu) * c2);
      }
    const have = ed.map(([i, j]) => hingeOf.get(key(i, j)) ?? -1);
    have.forEach((h, a) => (tri[3 * t + a] = h));
    const K = [0, 1, 2].filter((a) => have[a] >= 0);
    const B = [0, 1, 2].filter((a) => have[a] < 0);
    // condense: M' = M_KK - M_KB M_BB^-1 M_BK
    const inv = invertSmall(B.map((a) => B.map((b) => full[3 * a + b])));
    for (const a of K)
      for (const b of K) {
        let m = full[3 * a + b];
        for (let x = 0; x < B.length; x++)
          for (let y = 0; y < B.length; y++) m -= full[3 * a + B[x]] * inv[x][y] * full[3 * B[y] + b];
        M[9 * t + 3 * a + b] = m;
      }
  });
  return {
    w0: Int32Array.from(w0),
    w1: Int32Array.from(w1),
    e0: Int32Array.from(e0),
    e1: Int32Array.from(e1),
    sign,
    theta0,
    tri,
    M,
    sigma1: d.sigma1,
  };
}

function invertSmall(m: number[][]): number[][] {
  if (m.length === 0) return [];
  if (m.length === 1) return [[1 / m[0][0]]];
  if (m.length === 2) {
    const det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
    return [
      [m[1][1] / det, -m[0][1] / det],
      [-m[1][0] / det, m[0][0] / det],
    ];
  }
  // all three edges free: an isolated triangle has no bending energy
  return m.map((r) => r.map(() => 0));
}