// The whole precompute: from a cymsim configuration to the modal data the player needs.
//
// 1. profile, mesh (at the playback element size) and shell, as in the simulator;
// 2. static condensation of the in-plane and rotation DOFs (assemble.ts);
// 3. mass-normalised eigenmodes of the condensed transverse problem (eigen.ts);
// 4. per mode: the frequency, the quadratic form of the sigma1 loss, the complex listener gains
//    (gains.ts) and the mode shape on the nodes, quantised to int16.
//
// What is *not* baked in, so that the design panel can change it without a new eigensolve:
//
//   damping   sigma_m = sigma0 + sigma1 q_m with q_m = phi_m^T L phi_m stored as ``dampQ``.
//             q_m is invariant under a density change, because L is proportional to rho while a
//             mass-normalised phi is proportional to 1/sqrt(rho).
//   material  the condensed stiffness is proportional to E and the transverse mass to rho, both
//             exactly, so a change of either scales every frequency by sqrt((E/E0)/(rho/rho0))
//             and leaves the mode shapes unchanged apart from the normalisation factor
//             1/sqrt(rho/rho0). Poisson's ratio does change the shapes and needs a recompute.
//   listener  the gains are a sum over elements and can be redone in a fraction of a second.
//
// Everything else (radius, thickness, bell and bow geometry, imperfections, Poisson's ratio,
// mesh size) changes the eigenproblem and needs the full computation again.
import { cloneConfig } from "../config/yaml.ts";
import { defaultMallet, defaultStick, type SimulationConfig } from "../config/schema.ts";
import { CymbalProfile } from "../model/geometry.ts";
import { diskMesh } from "../model/mesh.ts";
import { buildShell } from "../model/shell.ts";
import { FORMAT_VERSION, type ModesData, type ModesMeta } from "../modes.ts";
import { condensedStiffness, type Progress } from "./assemble.ts";
import { symmetricEigen } from "./eigen.ts";
import { computeGains } from "./gains.ts";

export interface PrecomputeOptions {
  elementSize: number; // playback mesh element size [m]
  fmax: number; // keep modes below this frequency [Hz]
}

/** The parts of a configuration that the eigenproblem depends on. Two configurations with the
 * same structural key share a modes file; the rest is applied when the file is loaded. */
export function structuralKey(cfg: SimulationConfig, opts: PrecomputeOptions) {
  return {
    profile: cfg.profile,
    imperfections: cfg.imperfections,
    smoothing: cfg.mesh.smoothing,
    nu: cfg.material.nu,
    options: opts,
  };
}

export function computeModes(baseConfig: SimulationConfig, opts: PrecomputeOptions, progress: Progress): ModesData {
  const cfg = cloneConfig(baseConfig);
  cfg.mesh.element_size = opts.elementSize;
  progress("mesh and shell", 0);
  const profile = new CymbalProfile(cfg.profile, cfg.imperfections);
  const mesh = diskMesh(cfg.profile.radius, cfg.mesh.element_size, cfg.mesh.smoothing);
  const shell = buildShell(mesh, profile, cfg.material, cfg.numerics);
  const N = mesh.n, E = mesh.e, tri = mesh.tri;

  const S = condensedStiffness(mesh, shell, progress);

  // mass normalisation with the lumped transverse mass
  const Mw = new Float64Array(N), dinv = new Float64Array(N);
  for (let n = 0; n < N; n++) {
    Mw[n] = shell.M[5 * n + 2];
    dinv[n] = 1 / Math.sqrt(Mw[n]);
  }
  for (let i = 0; i < N; i++) for (let j = 0; j < N; j++) S[i * N + j] *= dinv[i] * dinv[j];
  progress("eigenmodes", 0);
  const eig = symmetricEigen(S, N, (phase, f) => progress(`eigenmodes: ${phase}`, f));

  // keep 5 Hz < f < fmax (drops the three rigid modes)
  const keep: number[] = [];
  for (let k = 0; k < N; k++) {
    const f = Math.sqrt(Math.max(eig.values[k], 0)) / (2 * Math.PI);
    if (f > 5 && f < opts.fmax) keep.push(k);
  }
  const Nm = keep.length;
  if (Nm === 0) throw new Error("no modes below the frequency limit: check the profile and the material");
  const freq = new Float64Array(Nm);
  const phi = new Float64Array(N * Nm); // node-major
  keep.forEach((k, m) => {
    freq[m] = Math.sqrt(eig.values[k]) / (2 * Math.PI);
    for (let n = 0; n < N; n++) phi[n * Nm + m] = eig.vectorsT[k * N + n] * dinv[n];
  });
  progress("damping and radiation", 0);

  // sigma1 quadratic form q_m = phi^T L phi with L = sum_e (rho h A)_e G^T G on w
  const dampQ = new Float64Array(Nm);
  const gxm = new Float64Array(Nm), gym = new Float64Array(Nm);
  for (let e = 0; e < E; e++) {
    gxm.fill(0);
    gym.fill(0);
    for (let i = 0; i < 3; i++) {
      const row = tri[3 * e + i] * Nm, gx = mesh.Gx[3 * e + i], gy = mesh.Gy[3 * e + i];
      for (let m = 0; m < Nm; m++) {
        gxm[m] += gx * phi[row + m];
        gym[m] += gy * phi[row + m];
      }
    }
    const c = shell.rhoHA[e];
    for (let m = 0; m < Nm; m++) dampQ[m] += c * (gxm[m] * gxm[m] + gym[m] * gym[m]);
  }

  // quantise the mode shapes
  const phiScale = new Float32Array(Nm), phiQ = new Int16Array(N * Nm);
  for (let m = 0; m < Nm; m++) {
    let top = 0;
    for (let n = 0; n < N; n++) top = Math.max(top, Math.abs(phi[n * Nm + m]));
    phiScale[m] = top / 32767 || 1;
  }
  for (let n = 0; n < N; n++) for (let m = 0; m < Nm; m++) phiQ[n * Nm + m] = Math.round(phi[n * Nm + m] / phiScale[m]);

  let mass = 0;
  for (const mw of Mw) mass += mw;
  const meta: ModesMeta = {
    format_version: FORMAT_VERSION,
    n_nodes: N, n_elements: E, n_modes: Nm, fmax: opts.fmax,
    radius: cfg.profile.radius, mass,
    E0: cfg.material.E, rho0: cfg.material.rho,
    config: baseConfig as unknown as Record<string, unknown>,
    playback_element_size: opts.elementSize,
    exciters: { stick: defaultStick(), mallet: defaultMallet() },
    notes: "phi = phi_q * phi_scale[m] is the transverse mode shape, mass-normalised at rho0, with " +
      "the in-plane and rotation DOFs statically condensed. Modal equation a'' + 2 sigma a' + " +
      "(2 pi f)^2 a = phi(p) . F, with sigma_m = sigma0 + sigma1 damp_q[m] and f scaled by " +
      "sqrt((E/E0)/(rho/rho0)). Pressure = Re(sum_m G_m * complex acceleration); unbaffled = front - back.",
  };

  const data: ModesData = {
    meta, nNodes: N, nElements: E, nModes: Nm, radius: cfg.profile.radius,
    xy: Float32Array.from(mesh.xy), tri: Int32Array.from(mesh.tri),
    z0: Float32Array.from(shell.z0Node), h: Float32Array.from(shell.hNode),
    freq, dampQ,
    gainFront: new Float32Array(2 * Nm), gainBack: new Float32Array(2 * Nm),
    phiScale, phiQ,
  };
  const g = computeGains(data, cfg.listener, freq, 1, undefined, (f) => progress("damping and radiation", f));
  data.gainFront = g.front;
  data.gainBack = g.back;
  progress("done", 1);
  return data;
}