/**
 * The computational grid: a cube, cell-centred, uniform in all three
 * directions.
 *
 * Fields are flattened x-fastest then y — the point (ix, iy, iz) is element
 * `ix + n*iy + n*n*iz` — which is the order the stencil shader indexes in and
 * the order the renderer reads a voxel in. Everything else treats a field as
 * an opaque npts-long array of f32.
 *
 * 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.
 *
 * `n` is required to be a multiple of 16 because the update shader dispatches
 * 8x8x4 workgroups and does not want a partial one at the far corner. That is
 * a convenience, not a physical constraint.
 */
export interface Grid {
  n: number;
  npts: number;
  /** Side length of the cubic domain, in metres, centred on the origin. */
  L: number;
  /** Grid spacing, L/n, in metres. */
  h: number;
}

export function makeGrid(n: number, L: number): Grid {
  return { n, npts: n * n * n, L, h: L / n };
}

/** The coordinate of cell index `i` along one axis, in metres. */
export const coord = (i: number, g: Grid): number => -g.L / 2 + (i + 0.5) * g.h;

/**
 * The timestep the explicit leapfrog is stable at, in 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) = 0.577. That is stricter than the 0.707 of the flat
 * case, which is the usual price of a dimension. `cfl` is the fraction of the
 * 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, cfl = 0.5): number {
  return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
}