import type { MethodDefinition } from "./types";
import {
  CURVED_EDGE_EQUATIONS,
  FREE_EDGE_EQUATIONS,
  GOVERNING_EQUATIONS,
  PINNED_CENTER_EQUATIONS,
  RADIATION_EQUATIONS,
  STRESS_EQUATIONS,
} from "./common";

export const IMPLICIT_METHOD: MethodDefinition = {
  id: "implicit",
  name: "Implicit finite differences",
  shortName: "Implicit FD",
  summary:
    "Same stencils as the explicit scheme, stiffness averaged in time; one sparse factorization, no step limit.",
  solverFunction: "plate_implicit",
  planFunction: "plate_implicit_plan",
  supports: {
    shape: ["rectangle", "ellipse"],
    bc: ["simply-supported", "clamped", "free", "pinned-center"],
    excitation: ["impulse", "harmonic"],
    audio: "clip",
  },
  readingOrder: [
    "plate_implicit.m",
    "plate_implicit_plan.m",
    "plate_domain.m",
    "plate_matrices.m",
    "plate_rigid_modes.m",
    "plate_curvature_free.m",
    "plate_curvature_domain.m",
    "plate_setup.m",
    "plate_force_time.m",
    "plate_force_footprint.m",
    "plate_stress.m",
    "plate_pad.m",
    "plate_curvature.m",
    "plate_listener.m",
    "plate_spectrum.m",
    "plate_result_schema.m",
  ],
  equations: [
    ...GOVERNING_EQUATIONS,
    {
      tex: String.raw`\frac{w^{n+1} - 2w^n + w^{n-1}}{\Delta t^2} + \sigma_0\frac{w^{n+1} - w^{n-1}}{\Delta t} - \sigma_1 \mathsf L\frac{w^{n+1} - w^{n-1}}{\Delta t} + \kappa^2 \mathsf B\left(\theta w^{n+1} + (1-2\theta) w^n + \theta w^{n-1}\right) = \frac{F^n s}{\mu}`,
      caption:
        "Theta scheme: the biharmonic term is averaged over three time levels. L and B are the sparse Laplacian and biharmonic matrices built from the same stencils as the explicit scheme, on the plate's footprint.",
    },
    {
      tex: String.raw`\mathbf A\, w^{n+1} = \mathbf C\, w^n - \mathbf A_{\mathrm{old}}\, w^{n-1} + \frac{\Delta t^2}{\mu} F^n s, \qquad \mathbf A = (1+\sigma_0\Delta t)\mathbf I - \sigma_1\Delta t\,\mathsf L + \theta\kappa^2\Delta t^2\,\mathsf B`,
      caption:
        "One linear system per step with a constant symmetric positive definite matrix: factor A once by sparse Cholesky, then two triangular solves per step.",
    },
    {
      tex: String.raw`\theta \ge \tfrac14 \;\Rightarrow\; \text{stable for every } \Delta t, \qquad \frac{T_{\text{num}}}{T} \approx 1 + \frac{(\omega\Delta t)^2}{12} \;\; (\theta = \tfrac14)`,
      caption:
        "No stability limit, so the step is chosen for accuracy. The price is period elongation: modes near half the sample rate run slow by tens of percent unless several steps per sample are used.",
    },
    ...FREE_EDGE_EQUATIONS,
    ...PINNED_CENTER_EQUATIONS,
    ...CURVED_EDGE_EQUATIONS,
    ...STRESS_EQUATIONS,
    ...RADIATION_EQUATIONS,
  ],
  validity: [
    "Second-order truncation error in dx and dy, identical to the explicit scheme, since the spatial stencils are the same.",
    "Second order in dt but dispersive: high modes are detuned rather than unstable. Two steps per sample keeps a 10 kHz mode within about 4 percent.",
    "The probe is snapped to the nearest grid node with a full stencil inside the plate; stress fields are evaluated only at the saved frames.",
    "Cost per step is a sparse triangular solve pair, about a millisecond for a 61 x 41 grid in the browser, independent of the plate stiffness.",
    "Audio covers only the simulated window.",
    "Free edges: the moment and shear conditions are natural conditions of a discrete bending energy, met to first order at the boundary and second order inside; the three rigid-body motions are removed after every step, so the displacement is that of a plate on soft supports, not one drifting in space.",
    "Pinned centre: the support is one grid node, so the stress right at it is resolution-dependent, and the two rigid rotations about it are removed rather than restrained by a modelled washer.",
    "Ellipse: the edge is a staircase of grid nodes, so the geometry, the edge conditions and the stress along the edge are first order in the grid spacing. The stress on the outermost ring of nodes is copied from the nearest node with a full stencil, which puts the peak stress of a clamped edge about one spacing inside. For simply supported edges the curvature term that the staircase loses is added back (the plate paradox), so the disk's frequencies converge to the exact ones.",
  ],
  durationOf: (spec) => spec.implicit.duration,
  gridOf: (spec) => [spec.implicit.nx, spec.implicit.ny],
  settingsKey: "implicit",
  usPerCellStep: 1.0,
};