Answer

The solution is complete and verified. Files created in ~: - ppm.js — the full implementation (609 lines), runs its own verification when executed directly - SOLUTION.md — the write-up below

Progressive Photon Mapping with a Median-Split kd-Tree, Cone Kernel, and Shrinking Radius

1. Problem and root-cause analysis

The class js-progressive-photon-mapping-kdtree-kernel-radius asks for a complete progressive photon mapper: photons emitted from lights, stored in a balanced kd-tree (median split on the widest axis), looked up with spherical / cone range queries, and used to estimate radiance with a cone-filtered kernel whose radius shrinks every pass according to the (n+1)/n accumulation rule. Surface and caustic maps must keep independent radii and accumulation state. The estimate must converge to the analytic diffuse-plane value, and caustic light behind a refractive sphere must be recovered from the kd-tree with the right photon counts and a monotonically decreasing radius.

The recurring root causes of failure in this class:

1. Linear scanning instead of a balanced kd-tree. A correct tree splits the current set at the median along the widest axis. Splitting on the first point or always on axis 0 destroys the O(log N) query and the balanced-depth guarantee.
2. Range/cone queries that drop boundary photons. Pruning must test d² <= R² on the split plane, and the node's own photon must be tested before descending. A cone query is a ball query (maxDist around the apex) plus the angular test dot(normalize(p-apex), axis) >= cos(halfAngle).
3. Wrong estimator normalization. Sampling a point light uniformly over the full sphere gives density n = P·cosθ/(4πd²) and photon power Φ_p = Φ/P, so n·Φ_p = Φ·cosθ/(4πd²) = E. The normalized flat cone kernel is K(d)=1/(πR²), hence L = (ρ/π)·(1/(πR²))·ΣΦ_p. Omitting the π in πR², or sampling only the downward hemisphere, silently doubles the density and yields a 2× error.
4. Radius that never shrinks. With constant R, extra passes only add noise. PPM accumulates N_{i+1} = N_i + M_i and keeps R²·N constant, so R_{i+1} = R_i·sqrt(N_i/N_{i+1}). For constant per-pass count this is exactly sqrt(n/(n+1)) = 1/sqrt((n+1)/n), i.e. R ∝ 1/sqrt(N).
5. Shared state between surface and caustic maps. Their densities differ, so each needs its own radius and accumulators.
6. Non-deterministic sampling. Use a seeded mulberry32 so the analytic check is reproducible.
7. Missing specular bookkeeping. Caustics are diffuse hits reached after a specular event. The tracer must refract with Snell's law, handle total internal reflection, and tag such hits.

2. Exact fix

Create ppm.js with the code below, then run node ppm.js.

'use strict';
/*
 * ppm.js - Progressive Photon Mapping (PPM) in JavaScript
 * =======================================================
 *
 * Features
 *   - Deterministic point-light / diffuse-plane / refractive-sphere scene
 *   - Balanced kd-tree photon map (median split on the widest axis)
 *   - Spherical range queries and apex/direction cone queries
 *   - Cone-filtered kernel radiance estimator
 *   - Progressive radius reduction:  N_{i+1} = N_i + alpha*M_i
 *                                    R_{i+1}^2 = R_i^2 * N_i / N_{i+1}
 *     (for a constant per-pass photon count M this is exactly the
 *      sqrt(n/(n+1)) = 1/sqrt((n+1)/n) accumulation rule)
 *   - Independent surface and caustic maps / accumulation state
 *   - Analytic diffuse-plane reference value
 *
 * Public API (CommonJS):
 *   const ppm = require('./ppm.js');
 *   ppm.PhotonMap, ppm.ProgressiveEstimator, ppm.Scene,
 *   ppm.PointLight, ppm.Plane, ppm.Sphere,
 *   ppm.tracePhoton, ppm.buildKdTree, ppm.rangeQuery, ppm.coneQuery,
 *   ppm.mulberry32, ppm.analyticPlaneRadiance, ppm.runVerification
 */

const PI = Math.PI;
const EPS = 1e-9;

/* ------------------------------------------------------------------ *
 * Deterministic PRNG (mulberry32)
 * ------------------------------------------------------------------ */
function mulberry32(seed) {
  let a = seed >>> 0;
  return function () {
    a = (a + 0x6d2b79f5) | 0;
    let t = Math.imul(a ^ (a >>> 15), 1 | a);
    t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

/* ------------------------------------------------------------------ *
 * Small vector helpers (plain arrays)
 * ------------------------------------------------------------------ */
const sub = (a, b) => [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
const add = (a, b) => [a[0] + b[0], a[1] + b[1], a[2] + b[2]];
const mul = (a, s) => [a[0] * s, a[1] * s, a[2] * s];
const dot = (a, b) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
const len = (a) => Math.sqrt(dot(a, a));
function normalize(a) {
  const l = len(a);
  return l > 0 ? [a[0] / l, a[1] / l, a[2] / l] : [0, 0, 0];
}
function dist2(a, b) {
  const dx = a[0] - b[0],
    dy = a[1] - b[1],
    dz = a[2] - b[2];
  return dx * dx + dy * dy + dz * dz;
}

/* ------------------------------------------------------------------ *
 * Balanced kd-tree (median split on the widest axis)
 * ------------------------------------------------------------------ */
function buildKdTree(photons, depth) {
  depth = depth || 0;
  if (!photons || photons.length === 0) return null;

  // Bounding box of the current subset.
  let mnx = Infinity, mny = Infinity, mnz = Infinity;
  let mxx = -Infinity, mxy = -Infinity, mxz = -Infinity;
  for (let i = 0; i < photons.length; i++) {
    const p = photons[i].position;
    if (p[0] < mnx) mnx = p[0];
    if (p[1] < mny) mny = p[1];
    if (p[2] < mnz) mnz = p[2];
    if (p[0] > mxx) mxx = p[0];
    if (p[1] > mxy) mxy = p[1];
    if (p[2] > mxz) mxz = p[2];
  }
  const wx = mxx - mnx, wy = mxy - mny, wz = mxz - mnz;
  let axis = 0, width = wx;
  if (wy > width) { axis = 1; width = wy; }
  if (wz > width) { axis = 2; width = wz; }

  const sorted = photons.slice().sort((a, b) => a.position[axis] - b.position[axis]);
  const mid = sorted.length >> 1;
  const node = {
    photon: sorted[mid],
    axis,
    depth,
    left: buildKdTree(sorted.slice(0, mid), depth + 1),
    right: buildKdTree(sorted.slice(mid + 1), depth + 1),
  };
  return node;
}

/* Spherical range query: all photons within `radius` of `center`. */
function rangeQuery(node, center, radius, out) {
  if (!node) return out;
  const p = node.photon.position;
  const dx = p[0] - center[0];
  const dy = p[1] - center[1];
  const dz = p[2] - center[2];
  if (dx * dx + dy * dy + dz * dz <= radius * radius) out.push(node.photon);

  const d = p[node.axis] - center[node.axis];
  if (d > 0) {
    rangeQuery(node.left, center, radius, out);
    if (d * d <= radius * radius) rangeQuery(node.right, center, radius, out);
  } else {
    rangeQuery(node.right, center, radius, out);
    if (d * d <= radius * radius) rangeQuery(node.left, center, radius, out);
  }
  return out;
}

/* Apex/direction cone query: photons inside `radius` of the apex whose
 * direction from the apex lies within `halfAngle` of `axis`. */
function coneQuery(node, apex, axis, halfAngle, radius, out) {
  if (!node) return out;
  const dir = normalize(axis);
  const cosHalf = Math.cos(halfAngle);

  // 1) kd-tree range query for the ball around the apex.
  const candidates = [];
  rangeQuery(node, apex, radius, candidates);

  // 2) angular filter.
  for (let i = 0; i < candidates.length; i++) {
    const ph = candidates[i];
    const w = sub(ph.position, apex);
    const l = len(w);
    if (l < EPS) { out.push(ph); continue; }
    const c = dot(w, dir) / l;
    if (c >= cosHalf) out.push(ph);
  }
  return out;
}

/* ------------------------------------------------------------------ *
 * Photon map wrapper
 * ------------------------------------------------------------------ */
class PhotonMap {
  constructor(photons) {
    this.photons = photons || [];
    this.root = buildKdTree(this.photons);
  }
  queryRadius(center, radius) {
    return rangeQuery(this.root, center, radius, []);
  }
  queryCone(apex, axis, halfAngle, radius) {
    return coneQuery(this.root, apex, axis, halfAngle, radius, []);
  }
  get size() {
    return this.photons.length;
  }
}

/* ------------------------------------------------------------------ *
 * Scene primitives
 * ------------------------------------------------------------------ */
class PointLight {
  constructor(position, power) {
    this.position = position;
    this.power = power;
  }
}

class Plane {
  constructor(y, albedo) {
    this.type = 'plane';
    this.y = y;
    this.albedo = albedo;
    this.normal = [0, 1, 0];
  }
  intersect(origin, dir) {
    if (Math.abs(dir[1]) < EPS) return null;
    const t = (this.y - origin[1]) / dir[1];
    return t > EPS ? t : null;
  }
}

class Sphere {
  constructor(center, radius, ior, albedo) {
    this.type = 'sphere';
    this.center = center;
    this.radius = radius;
    this.ior = ior;
    this.albedo = albedo == null ? 0 : albedo;
  }
  intersect(origin, dir) {
    const oc = sub(origin, this.center);
    const b = dot(oc, dir);
    const c = dot(oc, oc) - this.radius * this.radius;
    const disc = b * b - c;
    if (disc < 0) return null;
    const sq = Math.sqrt(disc);
    const t0 = -b - sq;
    const t1 = -b + sq;
    if (t0 > EPS) return t0;
    if (t1 > EPS) return t1;
    return null;
  }
}

class Scene {
  constructor(objects) {
    this.objects = objects;
  }
  /** Returns {t, object, point, normal} for the closest positive hit. */
  intersect(origin, dir) {
    let bestT = Infinity, bestObj = null;
    for (let i = 0; i < this.objects.length; i++) {
      const t = this.objects[i].intersect(origin, dir);
      if (t !== null && t < bestT) { bestT = t; bestObj = this.objects[i]; }
    }
    if (!bestObj) return null;
    const point = add(origin, mul(dir, bestT));
    let normal;
    if (bestObj.type === 'plane') normal = bestObj.normal;
    else normal = normalize(sub(point, bestObj.center));
    return { t: bestT, object: bestObj, point, normal };
  }
}

/* Refraction (Snell). `dir` and `normal` must be unit length.
 * Returns null on total internal reflection. */
function refract(dir, normal, eta) {
  let cosi = -dot(dir, normal);
  let n = normal;
  if (cosi < 0) { n = mul(normal, -1); cosi = -cosi; eta = 1 / eta; }
  const k = 1 - eta * eta * (1 - cosi * cosi);
  if (k < 0) return null;
  return normalize(add(mul(dir, eta), mul(n, eta * cosi - Math.sqrt(k))));
}

function reflect(dir, normal) {
  return normalize(sub(dir, mul(normal, 2 * dot(dir, normal))));
}

/* Uniform direction on the unit sphere (full solid angle). */
function sampleSphere(rng) {
  const z = 2 * rng() - 1;
  const phi = 2 * PI * rng();
  const s = Math.sqrt(Math.max(0, 1 - z * z));
  return [s * Math.cos(phi), z, s * Math.sin(phi)];
}

/**
 * Trace one photon through the scene.
 * Photons that reach a diffuse surface after >=1 specular (refractive)
 * interaction are tagged `caustic = true`.
 *
 * Returns the list of diffuse hits (position, power, incoming, caustic).
 */
function tracePhoton(scene, origin, dir, power, out, specularDepth) {
  out = out || [];
  specularDepth = specularDepth || 0;
  const hit = scene.intersect(origin, dir);
  if (!hit) return out;

  if (hit.object.type === 'plane') {
    out.push({
      position: hit.point,
      power,
      incoming: dir,
      normal: hit.object.normal,
      caustic: specularDepth > 0,
    });
    return out;
  }

  // Refractive sphere: transmit (or reflect on TIR), then continue.
  const n = hit.normal;
  const entering = dot(dir, n) < 0;
  const eta = entering ? 1 / hit.object.ior : hit.object.ior;
  let nd = refract(dir, entering ? n : mul(n, -1), eta);
  if (!nd) nd = reflect(dir, n);

  // Offset the origin slightly to avoid self-intersection.
  const off = add(hit.point, mul(nd, 1e-6));
  return tracePhoton(scene, off, nd, power, out, specularDepth + 1);
}

/**
 * Emit `count` photons from the light and trace them. Returns all diffuse
 * hits. Each photon carries power lightPower / count.
 */
function emitPhotons(scene, light, count, rng) {
  const power = light.power / count;
  const out = [];
  for (let i = 0; i < count; i++) {
    const dir = sampleSphere(rng);
    tracePhoton(scene, light.position, dir, power, out, 0);
  }
  return out;
}

/* ------------------------------------------------------------------ *
 * Progressive estimator
 * ------------------------------------------------------------------ */
class ProgressiveEstimator {
  /**
   * @param {object} opts
   *   position      hit point
   *   normal        surface normal (default [0,1,0])
   *   albedo        diffuse albedo
   *   radius        initial photon-gathering radius R0
   *   kernel        'flat' (indicator disc) | 'cone' (smooth cone filter)
   *   alpha         photon-accumulation factor (default 1)
   */
  constructor(opts) {
    this.position = opts.position;
    this.normal = opts.normal || [0, 1, 0];
    this.albedo = opts.albedo == null ? 0.6 : opts.albedo;
    this.initialRadius = opts.radius;
    this.radius = opts.radius;
    this.kernel = opts.kernel || 'flat';
    this.alpha = opts.alpha == null ? 1 : opts.alpha;

    this.passCount = 0;
    this.photonCount = 0; // N_i : accumulated photon count (radius rule)
    this.weight = 0;      // W   : sum of M_i used to average radiance
    this.flux = 0;        // S   : sum of M_i * I_i
    this.radiance = 0;    // current estimate
    this.history = [];
  }

  get brdf() {
    return this.albedo / PI;
  }

  /** Per-pass radiance estimate using the current radius. */
  estimatePass(photons) {
    const fr = this.brdf;
    const R2 = this.radius * this.radius;
    let acc = 0;
    if (this.kernel === 'cone') {
      for (let i = 0; i < photons.length; i++) {
        const d2 = dist2(photons[i].position, this.position);
        if (d2 > R2) continue;
        const K = (2 / (PI * R2)) * (1 - d2 / R2);
        acc += fr * photons[i].power * K;
      }
    } else {
      for (let i = 0; i < photons.length; i++) acc += fr * photons[i].power;
      acc /= PI * R2;
    }
    return acc;
  }

  /** Update with the photons found in this pass. Returns the new estimate. */
  update(photons) {
    this.passCount++;
    const M = photons.length;
    const I = this.estimatePass(photons);

    if (M > 0) {
      this.weight += M;
      this.flux += M * I;
    }
    if (this.weight > 0) this.radiance = this.flux / this.weight;

    // Radius accumulation rule: N_{i+1} = N_i + alpha*M_i and R^2 * N constant.
    const Nold = this.photonCount;
    const Nnew = Nold + this.alpha * M;
    if (Nold > 0) this.radius *= Math.sqrt(Nold / Nnew);
    this.photonCount = Nnew;

    this.history.push({
      pass: this.passCount,
      radius: this.radius,
      M,
      N: Nnew,
      radiance: this.radiance,
    });
    return this.radiance;
  }
}

/* ------------------------------------------------------------------ *
 * Analytic reference: diffuse plane under an isotropic point light
 *   L = rho * Phi * cos(theta) / (4 * pi^2 * d^2)
 * ------------------------------------------------------------------ */
function analyticPlaneRadiance(lightPos, point, normal, albedo, power) {
  const w = sub(lightPos, point);
  const d2 = dot(w, w);
  const d = Math.sqrt(d2);
  const cos = Math.max(0, dot(normalize(w), normal));
  return (albedo * power * cos) / (4 * PI * PI * d2);
}

/* ------------------------------------------------------------------ *
 * Verification suite
 * ------------------------------------------------------------------ */
function runVerification(opts) {
  opts = opts || {};
  const log = opts.log || console.log;
  const results = {};
  let ok = true;

  /* ---------------- Test 1: kd-tree correctness ---------------- */
  {
    const rng = mulberry32(7);
    const pts = [];
    for (let i = 0; i < 2000; i++) {
      pts.push({ position: [rng() * 4 - 2, rng() * 4 - 2, rng() * 4 - 2], power: 1 });
    }
    const map = new PhotonMap(pts);
    const center = [0.3, -0.2, 0.5];
    const radius = 0.7;
    const naive = pts.filter((p) => dist2(p.position, center) <= radius * radius);
    const q = map.queryRadius(center, radius);
    const countOk = q.length === naive.length;

    // Cone query vs brute force.
    const apex = [0, 0, 0];
    const axis = [1, 1, 0];
    const half = 0.5;
    const r2 = 1.2;
    const dir = normalize(axis);
    const cosHalf = Math.cos(half);
    const bruteCone = pts.filter((p) => {
      const w = sub(p.position, apex);
      const l = len(w);
      return l <= r2 && (l < EPS || dot(w, dir) / l >= cosHalf);
    });
    const cone = map.queryCone(apex, axis, half, r2);
    const coneOk = cone.length === bruteCone.length;

    // Balance: every leaf depth differs by at most 1 (median split guarantee).
    let minD = Infinity, maxD = -Infinity;
    (function walk(n) {
      if (!n) return;
      if (!n.left && !n.right) {
        minD = Math.min(minD, n.depth);
        maxD = Math.max(maxD, n.depth);
        return;
      }
      walk(n.left);
      walk(n.right);
    })(map.root);
    const balanced = maxD - minD <= 1;

    ok = ok && countOk && coneOk && balanced;
    results.kdTree = { rangeCount: q.length, naiveCount: naive.length, coneCount: cone.length, bruteCone: bruteCone.length, balanced, minDepth: minD, maxDepth: maxD };
    log(`[1] kd-tree: range ${q.length}/${naive.length}, cone ${cone.length}/${bruteCone.length}, balanced=${balanced} (depth ${minD}..${maxD})`);
  }

  /* ---------------- Test 2: diffuse plane convergence ---------------- */
  {
    const light = new PointLight([0, 2, 0], 100);
    const plane = new Plane(0, 0.6);
    const scene = new Scene([plane]);
    const hit = [1, 0, 0];
    const expected = analyticPlaneRadiance(light.position, hit, [0, 1, 0], 0.6, 100);

    const photonsPerPass = opts.photonsPerPass || 500000;
    const passes = opts.passes || 24;
    const radius0 = opts.radius0 || 0.1;

    const surfaceMapPhotons = [];
    const est = new ProgressiveEstimator({ position: hit, normal: [0, 1, 0], albedo: 0.6, radius: radius0, kernel: 'flat' });

    const seeds = opts.seeds || [12345];
    const perSeed = [];
    for (let s = 0; s < seeds.length; s++) {
      const rng = mulberry32(seeds[s]);
      const est2 = new ProgressiveEstimator({ position: hit, normal: [0, 1, 0], albedo: 0.6, radius: radius0, kernel: 'flat' });
      for (let p = 0; p < passes; p++) {
        const hits = emitPhotons(scene, light, photonsPerPass, rng).filter((h) => !h.caustic);
        const map = new PhotonMap(hits);
        const found = map.queryRadius(hit, est2.radius);
        const L = est2.update(found);
        surfaceMapPhotons.push(found.length);
      }
      perSeed.push({ L: est2.radiance, est: est2 });
    }
    const meanL = perSeed.reduce((a, b) => a + b.L, 0) / perSeed.length;
    const relErr = Math.abs(meanL - expected) / expected;
    const tol = opts.tolerance || 0.08;
    const convergeOk = relErr <= tol;
    const monotone = perSeed[0].est.history.every((h, i, arr) => i === 0 || h.radius <= arr[i - 1].radius + 1e-12);
    ok = ok && convergeOk && monotone;
    results.plane = {
      expected,
      estimate: meanL,
      relativeError: relErr,
      perSeed: perSeed.map((x) => x.L),
      monotoneRadius: monotone,
      firstRadius: perSeed[0].est.history[0].radius,
      lastRadius: perSeed[0].est.history[perSeed[0].est.history.length - 1].radius,
    };
    log(`[2] plane convergence: analytic=${expected.toFixed(5)} estimate=${meanL.toFixed(5)} relErr=${(relErr * 100).toFixed(2)}% (tol ${(tol * 100).toFixed(0)}%), monotoneRadius=${monotone}`);
  }

  /* ---------------- Test 3: caustic behind refractive sphere ---------------- */
  {
    const light = new PointLight([0, 3, 0], 100);
    const sphere = new Sphere([0, 0, 0], 1.0, 1.5, 0);
    const plane = new Plane(-2.5, 0.6);
    const scene = new Scene([sphere, plane]);

    const photonsPerPass = opts.causticPhotonsPerPass || 400000;
    const passes = opts.causticPasses || 12;
    const target = [0, -2.5, 0];
    const radius0 = opts.causticRadius0 || 0.35;

    const rng = mulberry32(4242);
    const est = new ProgressiveEstimator({ position: target, normal: [0, 1, 0], albedo: 0.6, radius: radius0, kernel: 'flat' });
    const counts = [];
    let causticTotal = 0;
    let lastMap = null;
    for (let p = 0; p < passes; p++) {
      const all = emitPhotons(scene, light, photonsPerPass, rng);
      const caustic = all.filter((h) => h.caustic);
      causticTotal += caustic.length;
      lastMap = new PhotonMap(caustic);
      const found = lastMap.queryRadius(target, est.radius);
      counts.push(found.length);
      est.update(found);
    }

    // Cone query from the light through the sphere toward the target.
    const apex = light.position;
    const axis = sub(target, apex);
    const coneHalf = 0.08;
    const coneDist = 6.5;
    const coneFound = lastMap.queryCone(apex, axis, coneHalf, coneDist);
    const dirCone = normalize(axis);
    const cosCone = Math.cos(coneHalf);
    const bruteCone = lastMap.photons.filter((ph) => {
      const w = sub(ph.position, apex);
      const l = len(w);
      return l <= coneDist && (l < EPS || dot(w, dirCone) / l >= cosCone);
    });
    const coneMatches = coneFound.length === bruteCone.length;
    const monotone = est.history.every((h, i, arr) => i === 0 || h.radius <= arr[i - 1].radius + 1e-12);
    const causticOk = causticTotal > 0 && counts[counts.length - 1] > 0 && monotone && coneMatches;

    ok = ok && causticOk;
    results.caustic = {
      totalCausticPhotons: causticTotal,
      perPassInRadius: counts,
      coneQueryCount: coneFound.length,
      bruteConeCount: bruteCone.length,
      coneMatchesBruteForce: coneMatches,
      finalRadiance: est.radiance,
      monotoneRadius: monotone,
      firstRadius: est.history[0].radius,
      lastRadius: est.history[est.history.length - 1].radius,
    };
    log(`[3] caustic: ${causticTotal} caustic photons, in-radius per pass=[${counts.join(',')}], coneQuery=${coneFound.length}/${bruteCone.length}, monotoneRadius=${monotone}, L=${est.radiance.toFixed(5)}`);
  }

  /* ---------------- Test 4: cone-filtered kernel convergence ---------------- */
  {
    const light = new PointLight([0, 2, 0], 100);
    const plane = new Plane(0, 0.6);
    const scene = new Scene([plane]);
    const hit = [1, 0, 0];
    const expected = analyticPlaneRadiance(light.position, hit, [0, 1, 0], 0.6, 100);
    const photonsPerPass = opts.photonsPerPass || 500000;
    const passes = opts.passes || 24;
    const rng = mulberry32(20240607);
    const est = new ProgressiveEstimator({ position: hit, normal: [0, 1, 0], albedo: 0.6, radius: opts.radius0 || 0.1, kernel: 'cone' });
    for (let p = 0; p < passes; p++) {
      const hits = emitPhotons(scene, light, photonsPerPass, rng).filter((h) => !h.caustic);
      const map = new PhotonMap(hits);
      const found = map.queryRadius(hit, est.radius);
      est.update(found);
    }
    const relErr = Math.abs(est.radiance - expected) / expected;
    const coneOk = relErr <= (opts.tolerance || 0.08) && est.history.every((h, i, arr) => i === 0 || h.radius <= arr[i - 1].radius + 1e-12);
    ok = ok && coneOk;
    results.coneKernel = { expected, estimate: est.radiance, relativeError: relErr };
    log(`[4] cone kernel: analytic=${expected.toFixed(5)} estimate=${est.radiance.toFixed(5)} relErr=${(relErr * 100).toFixed(2)}%`);
  }

  results.ok = ok;
  if (opts.log !== null) log(ok ? '\nALL CHECKS PASSED' : '\nSOME CHECKS FAILED');
  return results;
}

/* ------------------------------------------------------------------ */
module.exports = {
  PI,
  mulberry32,
  buildKdTree,
  rangeQuery,
  coneQuery,
  PhotonMap,
  PointLight,
  Plane,
  Sphere,
  Scene,
  refract,
  reflect,
  sampleSphere,
  tracePhoton,
  emitPhotons,
  ProgressiveEstimator,
  analyticPlaneRadiance,
  runVerification,
};

if (require.main === module) {
  runVerification();
}

3. Verification

$ node ppm.js
[1] kd-tree: range 52/52, cone 16/16, balanced=true (depth 10..10)
[2] plane convergence: analytic=0.27187 estimate=0.26998 relErr=0.69% (tol 8%), monotoneRadius=true
[3] caustic: 137344 caustic photons, in-radius per pass=[7194,7265,6310,6077,5899,5724,5470,5438,5311,5270,5252,5060], coneQuery=7888/7888, monotoneRadius=true, L=3.13828
[4] cone kernel: analytic=0.27187 estimate=0.28768 relErr=5.81%

ALL CHECKS PASSED

Multi-seed stress test:

$ node -e '
const ppm=require("./ppm.js");
const r=ppm.runVerification({seeds:[12345,999,7,42,2024,31337,1],passes:24,photonsPerPass:500000,log:()=>{}});
console.log("plane per-seed:",r.plane.perSeed.map(x=>x.toFixed(4)).join(" "));
console.log("expected",r.plane.expected.toFixed(4),"mean",r.plane.estimate.toFixed(4),
            "relerr",(r.plane.relativeError*100).toFixed(2)+"%");
console.log("caustic cone",r.caustic.coneQueryCount,"monotone",r.caustic.monotoneRadius,
            "radii",r.caustic.firstRadius.toFixed(4),"->",r.caustic.lastRadius.toFixed(4));
console.log("ok",r.ok);'
plane per-seed: 0.2700 0.2858 0.2740 0.2684 0.2690 0.2883 0.2778
expected 0.2719 mean 0.2762 relerr 1.58%
caustic cone 7888 monotone true radii 0.3500 -> 0.1120
ok true

Analytic reference: L = ρ·Φ·cosθ / (4π²d²) with ρ=0.6, Φ=100, light (0,2,0), hit (1,0,0) → L ≈ 0.271873.

Independence of the two maps. Test 2 builds the surface map from photons with caustic === false and its own estimator (R0 = 0.1). Test 3 builds the caustic map only from photons that refract through the glass sphere (caustic === true) with its own estimator (R0 = 0.35). Neither radius nor accumulator (N, W, S) is shared.