import { describe, expect, it } from "vitest";
import {
  FD_METHOD,
  SCENARIO_PRESETS,
  buildPlanBundle,
  cloneSpec,
  parseMethodPlan,
} from "../src/index";
import { maxAbs, runScript, runner, scalar, solve, testSpec, vector } from "./helpers";

describe("the hanging plate", () => {
  it("gets the compound-pendulum stiffness from the prestress operator to rounding", async () => {
    const spec = testSpec();
    const out = await runScript(
      spec,
      [
        "mdl = plate_model(params);",
        "d = norm([mdl.mnt.xc mdl.mnt.yc] - mdl.mnt.pivot);",
        "out.k_swing = mdl.R(:, 2)' * mdl.Kg * mdl.R(:, 2);",
        "out.mgd = mdl.mnt.mass * mdl.P.g * d;",
        "out.k_spin = mdl.R(:, 3)' * mdl.Kg * mdl.R(:, 3);",
        "out.k_heave = mdl.R(:, 1)' * mdl.Kg * mdl.R(:, 1);",
        "out.mass_check = mdl.ps.mass_check;",
        "out.kT = mdl.Krr(1, 1);",
        "out.T_over_L = mdl.mnt.mass * mdl.P.g / params.mount.length;",
        "out.hz = mdl.rigid_hz(:)';",
        "out.G = mdl.G(:)'; out.Krr = mdl.Krr(:)';",
      ].join("\n"),
    );
    const kSwing = scalar(out, "k_swing");
    const mgd = scalar(out, "mgd");
    expect(Math.abs(kSwing - mgd) / mgd).toBeLessThan(1e-10);
    expect(Math.abs(scalar(out, "k_spin"))).toBeLessThan(1e-9 * mgd);
    expect(Math.abs(scalar(out, "k_heave"))).toBeLessThan(1e-9 * mgd);
    expect(scalar(out, "mass_check")).toBeCloseTo(1, 10);
    expect(scalar(out, "kT")).toBeCloseTo(scalar(out, "T_over_L"), 8);
    // the two coupled pendulum modes of thread and plate, from the same 3 x 3 system
    const G = vector(out, "G");
    const K = vector(out, "Krr");
    const a = K[4] * G[0] - K[0] * G[4] * 0; // placeholder to keep the arrays used
    void a;
    const g11 = G[0],
      g12 = G[1],
      g22 = G[4];
    const k11 = K[0],
      k22 = K[4];
    // det(K - w^2 G) = 0 for the 2 x 2 sway-swing block
    const A = g11 * g22 - g12 * g12;
    const B = -(k11 * g22 + k22 * g11);
    const C = k11 * k22;
    const w2a = (-B - Math.sqrt(B * B - 4 * A * C)) / (2 * A);
    const w2b = (-B + Math.sqrt(B * B - 4 * A * C)) / (2 * A);
    const hz = Array.from(vector(out, "hz")).sort((p, q) => p - q);
    expect(hz[1]).toBeCloseTo(Math.sqrt(w2a) / (2 * Math.PI), 6);
    expect(hz[2]).toBeCloseTo(Math.sqrt(w2b) / (2 * Math.PI), 6);
    // a 300 x 200 mm steel sheet on a 300 mm thread: sway near 0.8 Hz, swing near 2.9 Hz
    expect(hz[1]).toBeGreaterThan(0.7);
    expect(hz[1]).toBeLessThan(0.9);
    expect(hz[2]).toBeGreaterThan(2.6);
    expect(hz[2]).toBeLessThan(3.2);
  });

  it("hangs with its centre of mass below an off-axis hole", async () => {
    const spec = testSpec();
    spec.mount = {
      type: "thread",
      x: 0.05,
      y: 0.19,
      length: 0.3,
      thread:
        spec.mount.type === "thread"
          ? spec.mount.thread
          : { name: "nylon", diameter: 0.5e-3, G: 1e9 },
    };
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(true);
    if (!run.ok) return;
    const plan = parseMethodPlan(run.value);
    // the hole snaps to the nearest grid node
    const hx = Math.round(0.05 / plan.dx) * plan.dx;
    const hy = Math.round(0.19 / plan.dy) * plan.dy;
    const dx = plan.xc - hx;
    const dy = plan.yc - hy;
    const n = Math.hypot(dx, dy);
    expect(plan.gx).toBeCloseTo(dx / n, 6);
    expect(plan.gy).toBeCloseTo(dy / n, 6);
    expect(plan.gz).toBe(0);
    expect(plan.tension.length).toBe(1);
    expect(plan.tension[0]).toBeCloseTo(plan.mass * 9.81, 6);
    expect(plan.sag).toBe(0);
  });

  it("is unstable when hung from below its centre of mass, and says so", async () => {
    const spec = testSpec();
    spec.mount = {
      type: "thread",
      x: 0.15,
      y: 0.01,
      length: 0.3,
      thread: { name: "nylon", diameter: 0.5e-3, G: 1e9 },
    };
    // the plate hangs from the hole in any case: plate_mount turns it so the mass is below the hole
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(true);
    if (!run.ok) return;
    const plan = parseMethodPlan(run.value);
    expect(plan.gy).toBeCloseTo(1, 6);
    expect(plan.rigid_hz.every((h) => Number.isFinite(h))).toBe(true);
  });

  it("shares the weight between two threads and restrains the yaw", async () => {
    const spec = testSpec();
    spec.mount = { type: "threads", x1: 0.06, y1: 0.19, x2: 0.24, y2: 0.19, length: 0.3 };
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(true);
    if (!run.ok) return;
    const plan = parseMethodPlan(run.value);
    expect(plan.tension.length).toBe(2);
    expect(plan.tension[0] + plan.tension[1]).toBeCloseTo(plan.mass * 9.81, 6);
    expect(plan.tension[0]).toBeCloseTo(plan.tension[1], 6);
    expect(plan.rigid_hz.every((h) => Number.isFinite(h) && h > 0.3)).toBe(true);
    expect(plan.rigid_names[2]).toMatch(/yaw/);
  });

  it("refuses two threads with the centre of mass outside them", async () => {
    const spec = testSpec();
    spec.mount = { type: "threads", x1: 0.2, y1: 0.19, x2: 0.28, y2: 0.19, length: 0.3 };
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(false);
    if (!run.ok) expect(run.error).toMatch(/slack/);
  });

  it("releases from a pull of the asked size and swings", async () => {
    const spec = testSpec();
    spec.excitation = { type: "release", x: 0.15, y: 0.02, displacement: 0.002, radius: 0.015 };
    spec.probe = { x: 0.15, y: 0.02 };
    spec.fd.duration = 0.001;
    const { result } = await solve(spec);
    expect(result.info.stable).toBe(true);
    // the probe node is near the pull point; the total departure there starts near the pull
    // (the probe sits a row above the pull's centre, where the swing is a little smaller)
    expect(Math.abs(result.probe.w[0])).toBeGreaterThan(0.0015);
    expect(Math.abs(result.probe.w[0])).toBeLessThan(0.0022);
    // the swing coordinate carries most of it and the spin none
    const c = result.rigid.series;
    expect(Math.abs(c[1][0])).toBeGreaterThan(0.005);
    expect(Math.abs(c[2][0])).toBeLessThan(1e-9);
    expect(result.info.notes.some((n) => /released/.test(n))).toBe(true);
  });
});

describe("plates on supports", () => {
  it("carries the weight on the pads at rest, sagging into them", async () => {
    const spec = testSpec("sheet-on-foam");
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(true);
    if (!run.ok) return;
    const plan = parseMethodPlan(run.value);
    expect(plan.gz).toBe(-1);
    expect(plan.sag).toBeGreaterThan(0);
    const { result } = await solve(spec);
    const total = result.supports.static.reduce((a, b) => a + b, 0);
    expect(total).toBeCloseTo(result.static.mass * 9.81, 6);
    expect(result.supports.static.every((f) => f > 0)).toBe(true);
    expect(result.static.sag).toBeGreaterThan(1e-4); // soft foam: a fraction of a millimetre
    expect(result.info.stable).toBe(true);
  });

  it("lifts off the pads under a hard strike and lands again", async () => {
    const spec = testSpec("sheet-on-foam");
    spec.excitation = {
      type: "impulse",
      x: 0.15,
      y: 0.1,
      force: 300,
      duration: 1e-3,
      radius: 0.02,
    };
    spec.fd.duration = 0.12;
    const { result } = await solve(spec);
    expect(result.info.stable).toBe(true);
    const f = result.supports.series[0];
    const n = f.length;
    const early = f.subarray(0, Math.floor(n / 3));
    const late = f.subarray(Math.floor((2 * n) / 3));
    expect(Math.min(...Array.from(early))).toBeLessThan(1e-6); // off the pad
    expect(Math.max(...Array.from(late))).toBeGreaterThan(result.supports.static[0] * 0.5); // back on it
  });

  it("holds the cymbal between two preloaded felts", async () => {
    const spec = testSpec("cymbal-on-stand");
    spec.fd.nx = 21;
    spec.fd.ny = 21;
    const { result } = await solve(spec);
    expect(result.supports.names).toEqual(["lower felt", "upper felt"]);
    const [lower, upper] = result.supports.static;
    expect(lower).toBeGreaterThan(0);
    expect(upper).toBeLessThan(0);
    expect(lower + upper).toBeCloseTo(result.static.mass * 9.81, 6);
    expect(result.info.stable).toBe(true);
  });

  it("sags like a cantilever beam in the vise", async () => {
    const spec = testSpec("cantilever-in-vise");
    spec.fd.nx = 31;
    spec.fd.ny = 7;
    spec.excitation = { type: "impulse", x: 0.2, y: 0.025, force: 1, duration: 1e-3, radius: 0.01 };
    const out = await runScript(
      spec,
      [
        "mdl = plate_model(params);",
        "w0 = reshape(mdl.st.w0, mdl.ny, mdl.nx);",
        "out.tip = w0(round((mdl.ny + 1) / 2), end);",
        "held = mdl.mnt.held; out.xclamp = max(mdl.dom.X(held));",
        "out.mu = mdl.P.mu; out.D = mdl.P.D; out.E = mdl.P.material.E; out.h = params.plate.h; out.Lx = params.plate.Lx;",
      ].join("\n"),
    );
    const Lfree = scalar(out, "Lx") - scalar(out, "xclamp");
    const q = scalar(out, "mu") * 9.81; // load per unit area
    const E = scalar(out, "E");
    const h = scalar(out, "h");
    // beam theory for a narrow strip: w = q L^4 / (8 E I / b) with I / b = h^3 / 12
    const beam = (q * Lfree ** 4) / (8 * E * (h ** 3 / 12));
    const tip = -scalar(out, "tip");
    expect(tip).toBeGreaterThan(0.85 * beam);
    expect(tip).toBeLessThan(1.15 * beam);
  });

  it("has no rigid-body coordinates in the vise and three elsewhere", async () => {
    const vise = await solve(testSpec("cantilever-in-vise"));
    expect(vise.result.rigid.nr).toBe(0);
    expect(vise.result.rigid.series.length).toBe(0);
    const hand = await solve(testSpec("disk-in-hand"));
    expect(hand.result.rigid.nr).toBe(3);
    expect(hand.result.rigid.hz.every((h) => Number.isFinite(h) && h >= 0)).toBe(true);
    expect(hand.result.rigid.hz.filter((h) => h > 0).length).toBeGreaterThanOrEqual(2);
  });
});

describe("the scheme", () => {
  it("keeps the swing of a struck hanging plate bounded and its split exact", async () => {
    const spec = testSpec();
    spec.fd.duration = 0.45;
    spec.damping.sigma0 = 0;
    spec.damping.sigma1 = 0;
    spec.output.nframes = 10;
    const out = await runScript(
      spec,
      [
        "S = plate_fd_setup(params); mdl = S.mdl;",
        "Rec = plate_fd_loop(S);",
        "E = zeros(3, S.nf); for k = 1:S.nf, E(:, k) = mdl.Q' * reshape(Rec.Wf(:, :, k), [], 1); end",
        "out.leak = max(max(abs(E(1:2, :))));",
        "out.swing = max(abs(Rec.C(2, :))); out.sway = max(abs(Rec.C(1, :)));",
        "n = size(Rec.C, 2); out.early = max(abs(Rec.C(2, 1:floor(n / 3)))); out.late = max(abs(Rec.C(2, floor(2 * n / 3):end)));",
        "% the angular impulse of the strike about the hole, and the swing it can give at most",
        "J = trapz(S.tt, S.F); rd = (params.excitation.x - mdl.mnt.pivot(1)) * mdl.mnt.ed(1) + (params.excitation.y - mdl.mnt.pivot(2)) * mdl.mnt.ed(2);",
        "out.bound = abs(J * rd) / mdl.G(2, 2) / (2 * pi * min(mdl.rigid_hz(mdl.rigid_hz > 0.1)));",
      ].join("\n"),
    );
    expect(scalar(out, "leak")).toBeLessThan(1e-12);
    // the swing stays of the order the impulse allows, and does not grow
    expect(scalar(out, "swing")).toBeLessThan(3 * scalar(out, "bound"));
    expect(scalar(out, "late")).toBeLessThan(
      1.5 * Math.max(scalar(out, "early"), scalar(out, "swing") * 0.5),
    );
  });

  it("factors every operator so the port can apply it as stencils and weights", async () => {
    const out = await runScript(
      testSpec(),
      [
        "S = plate_fd_setup(params);",
        "u = S.keep; Pz = spdiags(u, 0, S.N, S.N);",
        "KBf = Pz * (S.Gs' * S.WB * S.Gs) * Pz; out.eKB = full(max(max(abs(KBf - S.KB)))) / full(max(max(abs(S.KB))));",
        "KGf = Pz * (S.Gs' * S.WG * S.Gs) * Pz + spdiags(S.ks .* u, 0, S.N, S.N); out.eKG = full(max(max(abs(KGf - S.KG)))) / max(full(max(max(abs(S.KG)))), eps);",
        "Lf = spdiags(S.inva, 0, S.N, S.N) * (S.Gs' * S.WL * S.Gs) * Pz; out.eL = full(max(max(abs(Lf - S.Lm)))) / full(max(max(abs(S.Lm))));",
        "[gi, gj, gv] = find(S.Gs); out.eInt = max(abs(gv ./ S.Gscale(gi) - round(gv ./ S.Gscale(gi))));",
      ].join("\n"),
    );
    expect(scalar(out, "eKB")).toBeLessThan(1e-12);
    expect(scalar(out, "eKG")).toBeLessThan(1e-12);
    expect(scalar(out, "eL")).toBeLessThan(1e-12);
    expect(scalar(out, "eInt")).toBeLessThan(1e-9);
  });

  it("plans a time step at the stability bound", async () => {
    const spec = testSpec();
    const run = await runner.run({ bundle: buildPlanBundle(spec, FD_METHOD)! });
    expect(run.ok).toBe(true);
    if (!run.ok) return;
    const plan = parseMethodPlan(run.value);
    const { Lx, Ly, h } = spec.plate;
    const { E, nu, rho } = spec.material;
    const kappa = Math.sqrt((E * h ** 3) / (12 * (1 - nu ** 2)) / (rho * h));
    const dx = Lx / (spec.fd.nx - 1);
    const dy = Ly / (spec.fd.ny - 1);
    const lam = 4 / dx ** 2 + 4 / dy ** 2;
    const s1 = spec.damping.sigma1;
    const bending = 2 / (lam * (s1 + Math.sqrt(kappa ** 2 + s1 ** 2)));
    // the prestress and the thread spring lower the bound by a negligible amount
    expect(plan.dt_max).toBeLessThanOrEqual(bending);
    expect(plan.dt_max).toBeGreaterThan(0.999 * bending);
    expect(plan.dt).toBeCloseTo(spec.fd.cfl * plan.dt_max, 12);
    expect(plan.cell_steps).toBe(plan.ncells * plan.nsteps);
  });

  it("is stable just below the CFL limit and diverges just above it", async () => {
    const spec = testSpec();
    spec.fd.cfl = 0.98;
    const stable = await solve(spec);
    expect(stable.result.info.stable).toBe(true);
    expect(stable.result.probe.w.every(Number.isFinite)).toBe(true);

    spec.fd.cfl = 1.03;
    const blown = await solve(spec);
    expect(blown.result.info.stable).toBe(false);
    expect(blown.result.info.diverged_step).toBeGreaterThan(0);
    expect(blown.result.probe.w.every(Number.isFinite)).toBe(true);
    expect(blown.result.frames.w.data.every((v) => Number.isFinite(v) || Number.isNaN(v))).toBe(
      true,
    );
    expect(blown.result.info.notes.some((n) => /diverged/.test(n))).toBe(true);
  });

  it("runs every scenario on a coarse grid", async () => {
    for (const sc of SCENARIO_PRESETS) {
      const spec = testSpec(sc.id);
      const { result } = await solve(spec);
      expect(result.info.stable, sc.id).toBe(true);
      expect(maxAbs(result.probe.w), sc.id).toBeGreaterThan(0);
      expect(result.audio.p.every(Number.isFinite), sc.id).toBe(true);
      expect(result.frames.w.nf, sc.id).toBe(result.frames.t.length);
      expect(result.supports.series.length, sc.id).toBe(result.supports.names.length);
    }
  }, 240_000);

  it("reports its problem size and a series at the sample rate", async () => {
    const spec = testSpec();
    const { result } = await solve(spec);
    expect(result.info.ncells).toBe(spec.fd.nx * spec.fd.ny);
    expect(result.info.extra.cell_steps).toBe(result.info.ncells * result.info.nsteps);
    const expected = Math.floor(result.info.nsteps * result.info.dt * spec.output.fs);
    expect(result.t.length).toBe(expected);
    expect(result.audio.p.length).toBe(expected);
    expect(result.audio.model).toBe("unbaffled");
  });

  it("gives the same physics with the baffled listener except for the sound", async () => {
    const a = testSpec();
    const b = cloneSpec(a);
    b.listener.model = "baffled";
    const ra = await solve(a);
    const rb = await solve(b);
    expect(ra.result.probe.w).toEqual(rb.result.probe.w);
    expect(maxAbs(ra.result.audio.p)).toBeLessThan(maxAbs(rb.result.audio.p) * 1.5);
    expect(ra.result.audio.p).not.toEqual(rb.result.audio.p);
  });
});