/**
 * The computational grid: a square, cell-centred, uniform in both directions.
 *
 * Fields are flattened x-fastest — the point (ix, iy) is element `ix + nx*iy`
 * — which is the order the stencil shader indexes in and the order the
 * renderer reads a row of pixels in. Everything else in the project treats a
 * field as an opaque npts x 1 column vector.
 *
 * Cell-centred rather than node-centred so that no grid 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.
 */
export interface Grid {
  n: number;
  nx: number;
  ny: number;
  npts: number;
  /** Side length of the square domain, in metres, centred on the origin. */
  L: number;
  /** Grid spacing, L/n, in metres. */
  h: number;
  /** Coordinates of every point, npts each, x fastest — as the shaders see
   *  them (f32) and as the scene .m is evaluated at (f64). */
  x: Float32Array;
  y: Float32Array;
  x64: Float64Array;
  y64: Float64Array;
}

export function makeGrid(n: number, L: number): Grid {
  const h = L / n;
  const npts = n * n;
  const x64 = new Float64Array(npts);
  const y64 = new Float64Array(npts);
  for (let iy = 0; iy < n; iy++) {
    const yv = -L / 2 + (iy + 0.5) * h;
    for (let ix = 0; ix < n; ix++) {
      const k = ix + n * iy;
      x64[k] = -L / 2 + (ix + 0.5) * h;
      y64[k] = yv;
    }
  }
  return {
    n,
    nx: n,
    ny: n,
    npts,
    L,
    h,
    x: new Float32Array(x64),
    y: new Float32Array(y64),
    x64,
    y64,
  };
}

/**
 * The timestep the explicit leapfrog is stable at.
 *
 * All in SI: `h` in metres, `cmax` in metres per second, the result in
 * seconds.
 *
 * Leapfrog on p_tt = c^2 L p is stable while dt^2 c^2 |L|max <= 4, and the
 * discrete Laplacian's extreme eigenvalue is what differs between the
 * stencils: 8/h^2 for the 5-point one, 32/(3 h^2) for the 9-point
 * fourth-order one. That gives c*dt/h <= 1/sqrt(2) = 0.707 and
 * c*dt/h <= sqrt(3/8) = 0.612 respectively. `cfl` is the fraction of that
 * limit to run at, and `cmax` is the fastest sound speed anywhere in the
 * medium — a scatterer faster than the background sets the timestep for the
 * whole grid.
 */
export function stableDt(h: number, cmax: number, order: 2 | 4, cfl = 0.5): number {
  const limit = order === 2 ? Math.SQRT1_2 : Math.sqrt(3 / 8);
  return (cfl * limit * h) / Math.max(cmax, 1e-12);
}