/**
 * The two computational grids: a rectangular box of air, and a line of
 * string.
 *
 * The air grid is cell-centred and uniform, with the same spacing h along
 * every axis; the box is twice as long in x (the string direction) as in y
 * and z, so nx = 2*ny = 2*nz and a field is nx*ny*nz f32 values, flattened
 * x-fastest: the point (ix, iy, iz) is element `ix + nx*(iy + ny*iz)`.
 * Cell-centred so that no point sits exactly on the outer boundary: the
 * stencil takes the field outside the domain to be zero, and the absorbing
 * layer is meant to have swallowed the wave before it gets there.
 *
 * The string grid is node-centred: ns points from x = 0 to x = Ls inclusive,
 * spacing hs = Ls/(ns-1), so its endpoints land exactly on the terminations
 * that pin them. In the air's coordinates the string runs from -Ls/2 to
 * +Ls/2 along x, at y = 0 and a height above the top plate that the scene
 * decides.
 */
import { CFL } from './units.ts';

export interface AirGrid {
  nx: number;
  ny: number;
  nz: number;
  npts: number;
  /** Domain extents in metres, centred on the origin. */
  Lx: number;
  Ly: number;
  Lz: number;
  /** Grid spacing, the same in all three directions. */
  h: number;
  /** Coordinates of every grid point, npts each, x fastest — as the shaders
   *  see them (f32) and as the scene .m is evaluated at (f64). */
  x: Float32Array;
  y: Float32Array;
  z: Float32Array;
  x64: Float64Array;
  y64: Float64Array;
  z64: Float64Array;
}

export function makeAirGrid(nx: number, Lx: number): AirGrid {
  const ny = nx / 2;
  const nz = nx / 2;
  if (!Number.isInteger(ny)) throw new Error(`grid nx = ${nx} must be even`);
  const h = Lx / nx;
  const Ly = ny * h;
  const Lz = nz * h;
  const npts = nx * ny * nz;
  const x64 = new Float64Array(npts);
  const y64 = new Float64Array(npts);
  const z64 = new Float64Array(npts);
  for (let iz = 0; iz < nz; iz++) {
    const zv = -Lz / 2 + (iz + 0.5) * h;
    for (let iy = 0; iy < ny; iy++) {
      const yv = -Ly / 2 + (iy + 0.5) * h;
      const base = nx * (iy + ny * iz);
      for (let ix = 0; ix < nx; ix++) {
        const k = base + ix;
        x64[k] = -Lx / 2 + (ix + 0.5) * h;
        y64[k] = yv;
        z64[k] = zv;
      }
    }
  }
  return {
    nx, ny, nz, npts, Lx, Ly, Lz, h,
    x: new Float32Array(x64),
    y: new Float32Array(y64),
    z: new Float32Array(z64),
    x64, y64, z64,
  };
}

export interface StringGrid {
  ns: number;
  /** String length, metres. */
  Ls: number;
  /** Node spacing, Ls/(ns-1). */
  hs: number;
  /** Node positions, 0 to Ls inclusive. */
  xs: Float32Array;
  xs64: Float64Array;
  /** 1 at interior nodes, 0 at the two pinned ends. The compiled step
   *  multiplies its update by this, which is what terminates the string. */
  pin: Float32Array;
}

/**
 * The timestep the explicit leapfrog on the air grid is stable at, seconds.
 *
 * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
 * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
 * c*dt/h <= 1/sqrt(3). `cmax` is the fastest sound speed anywhere in the
 * medium; since the body's walls are a mask rather than a fast material (see
 * src/mgpu/ops.ts), cmax is normally just the speed of air.
 */
export function stableDt(h: number, cmax: number, cfl = CFL): number {
  return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
}

/**
 * The coarsest string spacing the stiff-string scheme demands at timestep dt,
 * in metres.
 *
 * The scheme (models/dulcimer.m) is the standard explicit leapfrog for
 *   u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
 * whose stability condition (Bilbao, Numerical Sound Synthesis, ch. 7) is
 *   hs^2 >= (a + sqrt(a^2 + 16 kap^2 dt^2)) / 2,  a = cs^2 dt^2 + 4 sig1 dt.
 */
export function stableHs(dt: number, cs: number, kap2: number, sig1: number): number {
  const a = cs * cs * dt * dt + 4 * sig1 * dt;
  return Math.sqrt((a + Math.sqrt(a * a + 16 * kap2 * dt * dt)) / 2);
}

/**
 * Build the string grid for length Ls at timestep dt, sized so that the
 * scheme is stable for every value the parameter sliders can reach.
 *
 * The air's CFL condition fixes dt, so the string has no say in the timestep;
 * what it gets to choose is its own spacing, and the finest stable spacing is
 * what maximizes the string's bandwidth. Sizing for the sliders' worst case
 * (highest fundamental, most stiffness, most frequency-dependent damping)
 * rather than their current values means moving a slider never forces a
 * recompile: the grid stays valid, merely a little coarser than that setting
 * alone would need.
 */
export function makeStringGrid(
  Ls: number,
  dt: number,
  worst: { f0: number; B: number; sig1: number },
): StringGrid {
  const cs = 2 * Ls * worst.f0;
  const kap2 = (worst.B * cs * cs * Ls * Ls) / (Math.PI * Math.PI);
  // 0.95: a little margin under the exact limit, as CFL is for the air.
  const hsMin = stableHs(dt, cs, kap2, worst.sig1) / 0.95;
  const ns = Math.max(8, Math.min(256, Math.floor(Ls / hsMin) + 1));
  const hs = Ls / (ns - 1);
  const xs64 = new Float64Array(ns);
  for (let i = 0; i < ns; i++) xs64[i] = i * hs;
  const pin = new Float32Array(ns).fill(1);
  pin[0] = 0;
  pin[ns - 1] = 0;
  return { ns, Ls, hs, xs: new Float32Array(xs64), xs64, pin };
}