/**
 * Scenes: what the medium is.
 *
 * A scene returns two fields of position, the sound speed c (m/s) and the
 * absorption sig (1/s). Both are ordinary fields, which is what lets one
 * function describe both the scatterer and the open boundary: the absorbing
 * layer around the outside is the statement that the medium swallows sound out
 * there, and a scene is free to put absorption inside the domain too, making a
 * lossy scatterer.
 *
 * These are built on the CPU, once, whenever a parameter changes. At 128^3
 * that is two million evaluations and takes well under a second; at 192^3 it
 * is seven million and is noticeable. The flat sibling writes its scenes as
 * MATLAB and evaluates them through numbl's interpreter, which is the more
 * interesting arrangement and the obvious thing to bring over here later.
 *
 * Interfaces are smoothed over about a cell. A jump between two neighbouring
 * cells is not resolved by the grid: it scatters the grid's own staircase
 * rather than the shape that was asked for.
 */
import { coord, type Grid } from './grid.ts';
import { C_AIR, SPONGE_FRAC, SPONGE_MAX } from './units.ts';

export interface SceneParam {
  key: string;
  label: string;
  min: number;
  max: number;
  step: number;
  value: number;
  unit?: string;
}

export interface Medium {
  c: Float32Array<ArrayBuffer>;
  sig: Float32Array<ArrayBuffer>;
  cmin: number;
  cmax: number;
}

export interface Scene {
  label: string;
  blurb: string;
  params: SceneParam[];
  /** Source settings this scene wants, applied when it is selected. */
  source?: Record<string, number>;
  /** Sound speed at (x, y, z), in m/s, and any absorption of its own. */
  medium(x: number, y: number, z: number, h: number, v: Record<string, number>): [number, number];
}

/**
 * Absorption profile for an open boundary: zero in the interior, ramping up
 * quadratically over a layer of width w inside each face of the cube and
 * reaching SPONGE_MAX at the wall. This is what makes the finite grid stand in
 * for an unbounded medium.
 *
 * The ramp is gradual on purpose. An absorbing layer is itself an impedance
 * mismatch, so a sudden one reflects; spreading it over a couple of
 * wavelengths keeps that small. It is not a perfectly matched layer, and at
 * grazing incidence it does leak.
 */
function sponge(x: number, y: number, z: number, L: number): number {
  const w = SPONGE_FRAC * L;
  const dx = Math.max(0, w - (L / 2 - Math.abs(x))) / w;
  const dy = Math.max(0, w - (L / 2 - Math.abs(y))) / w;
  const dz = Math.max(0, w - (L / 2 - Math.abs(z))) / w;
  const d = Math.max(dx, Math.max(dy, dz));
  return SPONGE_MAX * d * d;
}

/** A smoothed indicator: 1 well inside the surface, 0 well outside. */
const inside = (signedDistance: number, h: number): number =>
  0.5 * (1 - Math.tanh(signedDistance / (1.5 * h)));

/** Deterministic value noise on a coarse lattice, trilinearly interpolated. */
function makeNoise(m: number, seed: number): (u: number, v: number, w: number) => number {
  let s = (seed * 0x9e3779b9) >>> 0;
  const rand = () => {
    s = (s + 0x6d2b79f5) >>> 0;
    let t = s;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
  const g = new Float64Array(m * m * m);
  for (let i = 0; i < g.length; i++) g[i] = 2 * rand() - 1;
  const at = (i: number, j: number, k: number) => {
    const w = (a: number) => ((a % m) + m) % m;
    return g[w(i) + m * w(j) + m * m * w(k)];
  };
  // u, v, w in [0, 1) over the whole domain.
  return (u, v, w) => {
    const fx = u * m;
    const fy = v * m;
    const fz = w * m;
    const i = Math.floor(fx);
    const j = Math.floor(fy);
    const k = Math.floor(fz);
    const sx = fx - i;
    const sy = fy - j;
    const sz = fz - k;
    // Smoothstep on each axis so the field is continuous in its derivative,
    // which keeps the medium resolvable by the grid.
    const ex = sx * sx * (3 - 2 * sx);
    const ey = sy * sy * (3 - 2 * sy);
    const ez = sz * sz * (3 - 2 * sz);
    const lerp = (a: number, b: number, t: number) => a + (b - a) * t;
    const c00 = lerp(at(i, j, k), at(i + 1, j, k), ex);
    const c10 = lerp(at(i, j + 1, k), at(i + 1, j + 1, k), ex);
    const c01 = lerp(at(i, j, k + 1), at(i + 1, j, k + 1), ex);
    const c11 = lerp(at(i, j + 1, k + 1), at(i + 1, j + 1, k + 1), ex);
    return lerp(lerp(c00, c10, ey), lerp(c01, c11, ey), ez);
  };
}

/** Rebuilt whenever the random medium's seed or scale changes. */
let noiseCache: { key: string; f: (u: number, v: number, w: number) => number } | null = null;

export const scenes: Record<string, Scene> = {
  sphere: {
    label: 'Sphere',
    blurb:
      'One spherical scatterer in a uniform background: the reference case, and the ' +
      'three-dimensional problem with a classical series solution. Speeds above 1 behave ' +
      'nearly rigid, below 1 nearly pressure-release, and exactly 1 is no scatterer at all.',
    params: [
      { key: 'R', label: 'radius', min: 0.05, max: 0.5, step: 0.01, value: 0.25, unit: 'm' },
      { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 2.5 },
      { key: 'absorb', label: 'absorption', min: 0, max: 6000, step: 100, value: 0, unit: '1/s' },
    ],
    medium(x, y, z, h, v) {
      const r = Math.sqrt(x * x + y * y + z * z);
      const q = inside(r - v.R, h);
      return [C_AIR * (1 + (v.cin - 1) * q), v.absorb * q];
    },
  },

  pair: {
    label: 'Two spheres',
    blurb:
      'Two identical scatterers side by side. Each one reradiates what it receives, ' +
      'including what it receives from the other, so the pattern behind them is not the sum ' +
      'of two single-sphere patterns.',
    params: [
      { key: 'R', label: 'radius', min: 0.05, max: 0.35, step: 0.01, value: 0.16, unit: 'm' },
      { key: 'sep', label: 'separation', min: 0.15, max: 1.0, step: 0.01, value: 0.5, unit: 'm' },
      { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 3 },
    ],
    medium(x, y, z, h, v) {
      const d = v.sep / 2;
      const r1 = Math.sqrt(x * x + (y - d) * (y - d) + z * z);
      const r2 = Math.sqrt(x * x + (y + d) * (y + d) + z * z);
      const q = Math.max(inside(r1 - v.R, h), inside(r2 - v.R, h));
      return [C_AIR * (1 + (v.cin - 1) * q), 0];
    },
  },

  aperture: {
    label: 'Aperture',
    blurb:
      'A screen with a circular hole in it. Behind the hole is the three-dimensional ' +
      'diffraction pattern, which a plane cut through the flat problem cannot show: the ' +
      'Airy-like rings are a property of the circle, not of the slit.',
    params: [
      { key: 'a', label: 'hole radius', min: 0.04, max: 0.5, step: 0.01, value: 0.15, unit: 'm' },
      { key: 'th', label: 'screen thickness', min: 0.02, max: 0.2, step: 0.01, value: 0.06, unit: 'm' },
      { key: 'cw', label: 'screen speed', min: 0.05, max: 1, step: 0.05, value: 0.2 },
    ],
    medium(x, y, z, h, v) {
      // The screen is slow rather than fast. Reflection at an interface goes
      // as |c2 - c1|/(c2 + c1) at constant density, so c = 0.2 reflects about
      // as much as c = 5 would, but the timestep is set by the fastest speed
      // anywhere on the grid, so a slow screen is free and a fast one taxes
      // every step of the whole run. What a slow screen costs instead is
      // resolution inside itself, which its own absorption then swallows.
      const slab = inside(Math.abs(x) - v.th / 2, h);
      const rho = Math.sqrt(y * y + z * z);
      const hole = inside(rho - v.a, h);
      const q = slab * (1 - hole);
      return [C_AIR * (1 + (v.cw - 1) * q), 4000 * q];
    },
  },

  random: {
    label: 'Random medium',
    blurb:
      'Weak random structure everywhere, with no scatterer in particular. A pulse through ' +
      'it arrives on time and then keeps arriving: multiple scattering turns the tail into a ' +
      'coda, which is what a real inhomogeneous medium does.',
    params: [
      { key: 'amp', label: 'contrast', min: 0, max: 0.6, step: 0.01, value: 0.25 },
      { key: 'scale', label: 'blob size', min: 0.05, max: 0.6, step: 0.01, value: 0.18, unit: 'm' },
      { key: 'seed', label: 'seed', min: 1, max: 20, step: 1, value: 1 },
    ],
    medium(x, y, z, _h, v) {
      // `buildMedium` has rebuilt the lattice for this scale and seed before
      // calling us, so it is always here.
      const f = noiseCache!.f;
      const L = DOM.L;
      return [C_AIR * (1 + v.amp * f((x + L / 2) / L, (y + L / 2) / L, (z + L / 2) / L)), 0];
    },
  },
};

/** The domain the scenes are currently being evaluated over. Set by `build`;
 *  the random medium needs it to place its lattice. */
const DOM = { L: 2 };

export const sceneKeys = Object.keys(scenes);

/** Evaluate a scene over the grid, adding the absorbing layer. */
export function buildMedium(scene: Scene, grid: Grid, v: Record<string, number>): Medium {
  DOM.L = grid.L;
  if (scene === scenes.random) {
    const m = Math.max(2, Math.round(grid.L / Math.max(v.scale, 1e-3)));
    const key = `${m}:${v.seed}`;
    if (!noiseCache || noiseCache.key !== key) noiseCache = { key, f: makeNoise(m, v.seed) };
  }

  const { n, npts, L, h } = grid;
  const c = new Float32Array(npts);
  const sig = new Float32Array(npts);
  let cmin = Infinity;
  let cmax = -Infinity;
  for (let iz = 0; iz < n; iz++) {
    const z = coord(iz, grid);
    for (let iy = 0; iy < n; iy++) {
      const y = coord(iy, grid);
      const row = n * iy + n * n * iz;
      for (let ix = 0; ix < n; ix++) {
        const x = coord(ix, grid);
        const [cv, sv] = scene.medium(x, y, z, h, v);
        const k = ix + row;
        c[k] = cv;
        sig[k] = sv + sponge(x, y, z, L);
        if (cv < cmin) cmin = cv;
        if (cv > cmax) cmax = cv;
      }
    }
  }
  return { c, sig, cmin, cmax };
}