// Turning a selected linear mode into something the nonlinear simulator can run.
//
// Free decay: the shell is released from the mode shape at a chosen amplitude (given as the peak
// transverse displacement at the antinode, in units of the local thickness) and rings down. The
// in-plane displacements and rotations of the initial state come from the static condensation
// (the worker's "complement" request), so the released state is the complete linearised mode and
// not a pure-w field that would also launch fast membrane and rotation transients.
//
// Forcing: a distributed force F_n = Gamma * Mw_n * phi_n * env(t) * sin(omega_d t). Weighting by
// the transverse lumped mass makes the force M-orthogonal to every other mode, so in linear
// theory only the driven mode responds; whatever else appears in the nonlinear run got there
// through the coupling. Gamma is chosen so that the *linear* steady-state amplitude at resonance
// equals the requested one: Gamma = q_target * 2 sigma omega. With the damping scaled down this
// target is exceeded (at zero damping the linear response grows without bound), which is stated
// in the UI rather than hidden.
import type { ModesData } from "./modes.ts";
import type { DriveParams, ProjectionSpec } from "./model/simulator.ts";

export type Mechanism = "decay" | "forced";

export interface DriveSpec {
  mode: number;           // index into the modes file
  mechanism: Mechanism;
  amplitude: number;      // peak w at the antinode, in units of the local thickness
  detuneCents: number;    // forced: drive frequency offset from the mode frequency
  rampTime: number;       // forced: seconds of linear onset
  driveTime: number;      // forced: drive stops after this many seconds (free decay afterwards)
}

/** Mode shape m on the nodes, dequantised, times ``factor`` (the density rescaling phiFactor). */
export function modeShape(modes: ModesData, m: number, factor = 1): Float64Array {
  const N = modes.nNodes, Nm = modes.nModes, out = new Float64Array(N);
  const s = modes.phiScale[m] * factor;
  for (let n = 0; n < N; n++) out[n] = modes.phiQ[n * Nm + m] * s;
  return out;
}

/** Node of largest |phi|. */
export function antinode(shape: Float64Array): { node: number; absPhi: number } {
  let node = 0, absPhi = 0;
  for (let n = 0; n < shape.length; n++) if (Math.abs(shape[n]) > absPhi) { absPhi = Math.abs(shape[n]); node = n; }
  return { node, absPhi };
}

/** Modal-coordinate amplitude that puts the antinode at ``woverh`` thicknesses. */
export function modalAmplitude(modes: ModesData, shape: Float64Array, woverh: number): number {
  const { node, absPhi } = antinode(shape);
  return (woverh * modes.h[node]) / absPhi;
}

/** Initial state (q, v) of the free decay: w = a phi, in-plane and rotations from the complement
 * o (node-major u, v, thx, thy per node, computed for the *unit* shape), velocities zero. */
export function modalInitialState(shape: Float64Array, o: Float64Array, a: number): { q: Float64Array; v: Float64Array } {
  const N = shape.length, q = new Float64Array(5 * N), v = new Float64Array(5 * N);
  for (let n = 0; n < N; n++) {
    q[5 * n] = a * o[4 * n];
    q[5 * n + 1] = a * o[4 * n + 1];
    q[5 * n + 2] = a * shape[n];
    q[5 * n + 3] = a * o[4 * n + 2];
    q[5 * n + 4] = a * o[4 * n + 3];
  }
  return { q, v };
}

/** The distributed drive of the forced mechanism. ``Mw`` is the transverse lumped mass of the
 * simulator's shell (current material), ``sigma`` the decay rate used for the amplitude target,
 * ``omegaMode`` the (current-material) angular frequency of the mode. */
export function buildDrive(spec: DriveSpec, modes: ModesData, shape: Float64Array, Mw: Float64Array,
                           sigma: number, omegaMode: number, dt: number): DriveParams {
  const qTarget = modalAmplitude(modes, shape, spec.amplitude);
  const gamma = qTarget * 2 * Math.max(sigma, 1e-6) * omegaMode;
  const N = shape.length, fNode = new Float64Array(N);
  for (let n = 0; n < N; n++) fNode[n] = gamma * Mw[n] * shape[n];
  const omega = omegaMode * Math.pow(2, spec.detuneCents / 1200);
  return {
    fNode, omega,
    rampSteps: Math.max(1, Math.round(spec.rampTime / dt)),
    stopStep: Math.max(1, Math.round(spec.driveTime / dt)),
  };
}

/** The projection matrix (mode-major, mass-weighted): with mass-normalised shapes and the same
 * shell mass, q_m = sum_n phiM[m N + n] w_n is the modal coordinate. ``Mw`` must be the
 * transverse lumped mass of the shell actually being simulated. */
export function buildProjection(modes: ModesData, Mw: Float64Array, factor: number, every: number): ProjectionSpec {
  const N = modes.nNodes, Nm = modes.nModes;
  const phiM = new Float32Array(Nm * N);
  for (let m = 0; m < Nm; m++) {
    const s = modes.phiScale[m] * factor, row = m * N;
    for (let n = 0; n < N; n++) phiM[row + n] = modes.phiQ[n * Nm + m] * s * Mw[n];
  }
  return { phiM, nModes: Nm, every };
}