// Contact patch of one hit, resolved on the modal data: patch points (a line for a stick, a
// centre plus hexagon for a mallet), their containing triangles, the residual heights above the
// exciter's best-fit face plane, and the mode shape values at the points. Mirrors
// cymsim/exciter.py build_patch, with the surface height taken from the nodal z0 by linear
// interpolation.
import type { Exciter } from "./config/schema.ts";
import { planeResiduals } from "./model/linalg.ts";
import { interpolateNodal, locate, modeValuesAt, type ModesData } from "./modes.ts";

export interface HitPatch {
  phi: Float32Array; // (P, Nm) point-major mode values
  offsets: Float32Array; // (P,)
  kPoint: number;
  alpha: number;
  mass: number;
  gap: number;
}

export function contactStiffness(ex: Exciter, velocity: number): number {
  if (ex.kind === "stick") return ex.k;
  return ((ex.alpha + 1) * 0.5 * ex.mass * velocity * velocity) / Math.pow(ex.compression, ex.alpha + 1);
}

/** Half-extent of the patch [m], used to keep hits inside the mesh. */
export function patchExtent(ex: Exciter): number {
  return ex.kind === "stick" ? ex.patch_length / 2 : ex.radius;
}

/** ``maxPoints`` limits the number of springs used in real time (the exported configuration keeps
 * the exciter's own patch): a stick keeps the centre and the ends of its line, a mallet its centre. */
export function buildHitPatch(data: ModesData, ex: Exciter, r: number, phiAngle: number, velocity: number,
                              maxPoints = 7, phiFactor = 1): HitPatch | null {
  const xs = r * Math.cos(phiAngle), ys = r * Math.sin(phiAngle);
  const px: number[] = [], py: number[] = [];
  if (ex.kind === "stick") {
    const ang = phiAngle + ex.patch_angle, P = Math.max(1, Math.min(Math.floor(ex.patch_points), maxPoints));
    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));
    }
  } else {
    px.push(xs);
    py.push(ys);
    for (let k = 0; k < 6 && maxPoints > 1; 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, Nm = data.nModes;
  const phi = new Float32Array(P * Nm), z0p = new Float64Array(P);
  for (let p = 0; p < P; p++) {
    const loc = locate(data, px[p], py[p]);
    if (!loc) return null;
    modeValuesAt(data, loc.nodes, loc.weights, phi, p * Nm, phiFactor);
    z0p[p] = interpolateNodal(data.z0, loc.nodes, loc.weights);
  }
  const offsets = Float32Array.from(planeResiduals(px.map((x) => x - xs), py.map((y) => y - ys), z0p));
  return { phi, offsets, kPoint: contactStiffness(ex, velocity) / P, alpha: ex.alpha, mass: ex.mass, gap: ex.gap };
}