The spreading laboratory

export function gaussianPropagator(x: number, t: number, D: number): Evaluation {
  if (!Number.isFinite(x) || !validDt(D, t))
    return outside(
      "probabilityDensity",
      "gaussianPropagator",
      "finite x, D >= 0, t >= 0",
      "Use finite position and nonnegative diffusivity and time.",
    );
  if (D === 0 || t === 0) return pointMass("probabilityDensity", "gaussianPropagator");
  const sigma = scale(D, t);
  const z = x / sigma;
  if (!(sigma > 0) || !Number.isFinite(sigma))
    return outside(
      "probabilityDensity",
      "gaussianPropagator",
      "representable spread",
      "The spread is outside the numerical range.",
      "numerical",
    );
  return number(
    "probabilityDensity",
    "gaussianPropagator",
    Math.exp(-Math.log(sigma) - 0.5 * Math.log(2 * Math.PI) - 0.5 * z * z),
    undefined,
    undefined,
    true,
  );
}

/** Closed intervals include the atom at t=0, including the interval [0,0]. */
export function intervalProbability(x1: number, x2: number, t: number, D: number): Evaluation {
  const owner = "intervalProbability";
  if (![x1, x2].every(Number.isFinite) || x1 > x2 || !validDt(D, t))
    return outside(
      "intervalProbability",
      owner,
      "x1 <= x2, D >= 0, t >= 0",
      "Use ordered finite interval endpoints and nonnegative diffusivity and time.",
    );
  if (D === 0 || t === 0) return number("intervalProbability", owner, x1 <= 0 && x2 >= 0 ? 1 : 0);
  if (x1 === x2) return number("intervalProbability", owner, 0);
  const s = 2 * Math.sqrt(D) * Math.sqrt(t);
  if (!(s > 0) || !Number.isFinite(s))
    return outside(
      "intervalProbability",
      owner,
      "representable spread",
      "The spread is outside the numerical range.",
      "numerical",
    );
  const lo = x1 / s;
  const hi = x2 / s;
  const width = (x2 - x1) / s;
  const mid = lo + width / 2;
  let p: number;
  // Very narrow intervals need direct integration, even when both tails are small.
  if (width > 0 && width * (1 + Math.abs(mid)) < 0.01) {
    const nodes = [0.1834346424956498, 0.525532409916329, 0.7966664774136267, 0.9602898564975363];
    const weights = [0.362683783378362, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763];
    let sum = 0;
    for (let i = 0; i < 4; i++) {
      const offset = (width / 2) * nodes[i]!;
      sum += weights[i]! * (Math.exp(-((mid - offset) ** 2)) + Math.exp(-((mid + offset) ** 2)));
    }
    p = (width / (2 * Math.sqrt(Math.PI))) * sum;
  } else
    p =
      lo >= 0
        ? 0.5 * (tail(lo) - tail(hi))
        : hi <= 0
          ? 0.5 * (tail(-hi) - tail(-lo))
          : 0.5 * (central(hi) - central(lo));
  if (p < 0 || p > 1)
    return outside(
      "intervalProbability",
      owner,
      "probability range",
      "Numerical evaluation did not produce an admissible probability.",
      "numerical",
    );
  return number("intervalProbability", owner, p, undefined, undefined, true);
}

export function erf(x: number): number {
  validate(x);
  if (x < 0) return -erf(-x);
  return x < 1.5 ? erfPositive(x) : 1 - erfcPositive(x);
}

export function erfc(x: number): number {
  validate(x);
  return x < 0 ? 2 - erfcPositive(-x) : erfcPositive(x);
}

export function radialPropagator2d(r: number, t: number, D: number): Evaluation {
  return radial(r, t, D, 2);
}

export function radialPropagator3d(r: number, t: number, D: number): Evaluation {
  return radial(r, t, D, 3);
}

export function ftcs1d(p: FtcsParameters, budget = REFERENCE_FTCS_BUDGET): Computation<FtcsFrames> {
  const { n, frames, stepsPerFrame, D, dx, dt, profile } = p;
  if (!integer(n, 3) || !integer(frames, 1) || !integer(stepsPerFrame, 1))
    return invalid(
      ["n", "frames", "stepsPerFrame"],
      "Use at least three cells and positive whole-number frame and step counts.",
    );
  if (![D, dx, dt].every(Number.isFinite))
    return {
      kind: "refused",
      refusal: makeRefusal("nonfinite-input", { parameterIds: ["D", "dx", "dt"] }),
    };
  if (D < 0 || dx <= 0 || dt <= 0 || ![0, 1, 2].includes(profile))
    return invalid(
      ["D", "dx", "dt", "profile"],
      "Use nonnegative D, positive dx and dt, and profile 0, 1, or 2.",
    );
  // This exact operation order is part of the planned upstream conformance contract.
  const r = D === 0 ? 0 : (D * dt) / (dx * dx);
  if (!Number.isFinite(r) || (D > 0 && r === 0))
    return numerical("The stability ratio is outside binary64 range; the grid was not advanced.");
  if (r > 0.5) {
    const dtMax = (dx * dx) / (2 * D);
    const dxMin = Math.sqrt(2 * D * dt);
    const dMax = (dx * dx) / (2 * dt);
    if (![dtMax, dxMin, dMax].every((v) => Number.isFinite(v) && v > 0))
      return numerical("The dimensional stability repairs are outside binary64 range.");
    // Mathematical bounds can round to a value with r=0.5000000000000001.
    // Preserve dtMax in the explanation; make each offered action executable.
    const dtRepair = stableRepair(dtMax, false, (v) => (D * v) / (dx * dx));
    const dxRepair = stableRepair(dxMin, true, (v) => (D * dt) / (v * v));
    const dRepair = stableRepair(dMax, false, (v) => (v * dt) / (dx * dx));
    if (dtRepair === null || dxRepair === null || dRepair === null)
      return numerical(
        "No representable repair was found near the dimensional stability boundary.",
      );
    return {
      kind: "refused",
      refusal: makeRefusal(
        "ftcs-unstable",
        { parameterIds: ["dt", "dx", "D"] },
        {
          details: { ratio: r, limit: 0.5, dtMax },
          rankedRepairs: [
            {
              label: "Reduce the time step to the explicit scheme's limit.",
              action: { parameterId: "dt", value: dtRepair },
            },
            {
              label: "Use a coarser spatial grid.",
              action: { parameterId: "dx", value: dxRepair },
            },
            {
              label: "Choose a smaller diffusivity; this changes the physical setup.",
              action: { parameterId: "D", value: dRepair },
            },
          ],
        },
      ),
    };
  }
  const stepCount = (frames - 1) * stepsPerFrame;
  const limited = budgetCheck(n * stepCount, 8 * n * (frames + 2), budget);
  if (limited) return limited;
  const elapsedTime = stepCount * dt;
  if (!Number.isFinite(elapsedTime) || (stepCount > 0 && elapsedTime === 0))
    return numerical("Elapsed model time is outside binary64 range.");
  const field = new Float64Array(n);
  if (profile === 0) field[Math.floor(n / 2)] = 1 / dx;
  else if (profile === 1) field.fill(1, 0, Math.floor(n / 2));
  else {
    field[Math.floor(n / 4)] = 0.5 / dx;
    field[Math.floor((3 * n) / 4)] = 0.5 / dx;
  }
  if (!field.every(Number.isFinite))
    return numerical("The initial cell density is outside binary64 range.");
  const derivative = new Float64Array(n);
  const values = new Float64Array(n * frames);
  values.set(field);
  for (let frame = 1; frame < frames; frame++) {
    if (!advanceInPlace(field, derivative, r, stepsPerFrame))
      return numerical("A cell became nonfinite or negative; no partial frames are returned.");
    values.set(field, frame * n);
  }
  return {
    kind: "accepted",
    data: Object.freeze({
      values,
      shape: Object.freeze([frames, n] as const),
      stepCount,
      elapsedTime,
      stabilityRatio: r,
      boundary: "zero-flux",
      ownerId: "diffusion.ftcs1d",
    }),
  };
}