// Exciters and contact patches, as in cymsim/exciter.py: a rigid mass on P power-law springs
// F = k delta^alpha along a line (stick) or on a small disk (mallet), released at its hit time
// just above the moving surface. All hits are stacked into arrays padded to the largest patch.
import type { Hit } from "../config/schema.ts";
import type { CymbalProfile } from "./geometry.ts";
import { planeResiduals } from "./linalg.ts";
import { locate, type Mesh } from "./mesh.ts";

export const HOLD_HEIGHT = 1e3;

export interface ContactPatch {
  hit: Hit; x: Float64Array; y: Float64Array; nodes: Uint32Array; weights: Float64Array; offsets: Float64Array;
  kTotal: number; alpha: number; mass: number; gap: number; velocity: number; time: number;
}

export function patchEnergy(p: ContactPatch): number { return 0.5 * p.mass * p.velocity * p.velocity; }
export function maxCompression(p: ContactPatch): number { return Math.pow(((p.alpha + 1) * patchEnergy(p)) / p.kTotal, 1 / (p.alpha + 1)); }

/** Contact time against a rigid wall, 2 c(alpha) d_max / v with c = int_0^1 dx / sqrt(1 - x^(alpha+1)). */
export function freeContactTime(p: ContactPatch): number {
  // substitute x = 1 - s^2 to remove the endpoint singularity, then Simpson on s in [0, 1]
  const a1 = p.alpha + 1, n = 2000;
  const f = (s: number) => (s === 0 ? 2 / Math.sqrt(a1) : (2 * s) / Math.sqrt(1 - Math.pow(1 - s * s, a1)));
  let sum = f(0) + f(1);
  for (let i = 1; i < n; i++) sum += (i % 2 ? 4 : 2) * f(i / n);
  const c = sum / (3 * n);
  return (2 * c * maxCompression(p)) / p.velocity;
}

export function buildPatch(hit: Hit, mesh: Mesh, profile: CymbalProfile): ContactPatch {
  const ex = hit.exciter;
  const xs = hit.r * Math.cos(hit.phi), ys = hit.r * Math.sin(hit.phi);
  let px: number[], py: number[], kTotal: number, alpha: number;
  if (ex.kind === "stick") {
    const ang = hit.phi + ex.patch_angle, P = Math.max(1, Math.floor(ex.patch_points));
    px = []; py = [];
    for (let i = 0; i < P; i++) {
      const xi = (P === 1 ? 0 : -0.5 + i / (P - 1)) * ex.patch_length;
      px.push(xs + xi * Math.cos(ang)); py.push(ys + xi * Math.sin(ang));
    }
    kTotal = ex.k; alpha = ex.alpha;
  } else {
    alpha = ex.alpha;
    kTotal = ((alpha + 1) * 0.5 * ex.mass * hit.velocity ** 2) / ex.compression ** (alpha + 1);
    px = [xs]; py = [ys];
    for (let k = 0; k < 6; k++) { const th = (k * Math.PI) / 3; px.push(xs + ex.radius * Math.cos(th)); py.push(ys + ex.radius * Math.sin(th)); }
  }
  const P = px.length, nodes = new Uint32Array(3 * P), weights = new Float64Array(3 * P), z0p = new Float64Array(P);
  for (let p = 0; p < P; p++) {
    const loc = locate(mesh, px[p], py[p]);
    nodes.set(loc.nodes, 3 * p); weights.set(loc.weights, 3 * p);
    z0p[p] = profile.z0(px[p], py[p]);
  }
  const offsets = planeResiduals(px.map((x) => x - xs), py.map((y) => y - ys), z0p);
  return { hit, x: Float64Array.from(px), y: Float64Array.from(py), nodes, weights, offsets, kTotal, alpha, mass: ex.mass,
           gap: ex.gap, velocity: hit.velocity, time: hit.time };
}

export interface HitSchedule {
  patches: ContactPatch[]; K: number; P: number;
  nodes: Uint32Array; weights: Float64Array; offsets: Float64Array;   // (K, P, 3), (K, P, 3), (K, P)
  kPoint: Float64Array; alpha: Float64Array; mass: Float64Array; gap: Float64Array; velocity: Float64Array;
  releaseStep: Uint32Array;
}

export function buildSchedule(patches: ContactPatch[], dt: number): HitSchedule {
  // K = 0 (no exciters at all) is a valid schedule here: cymmodes drives the shell through
  // initial conditions and distributed forces instead of contacts.
  const K = patches.length, P = K ? Math.max(...patches.map((p) => p.x.length)) : 0;
  const nodes = new Uint32Array(K * P * 3), weights = new Float64Array(K * P * 3), offsets = new Float64Array(K * P).fill(-1e6);
  patches.forEach((p, k) => {
    const n = p.x.length;
    nodes.set(p.nodes, (k * P) * 3); weights.set(p.weights, (k * P) * 3); offsets.set(p.offsets, k * P);
    void n;
  });
  return {
    patches, K, P, nodes, weights, offsets,
    kPoint: Float64Array.from(patches, (p) => p.kTotal / p.x.length),
    alpha: Float64Array.from(patches, (p) => p.alpha),
    mass: Float64Array.from(patches, (p) => p.mass),
    gap: Float64Array.from(patches, (p) => p.gap),
    velocity: Float64Array.from(patches, (p) => p.velocity),
    releaseStep: Uint32Array.from(patches, (p) => Math.round(p.time / dt)),
  };
}