/**
 * The density field.
 *
 * Turns the projected manifold — 34 points on a plane — into a continuous
 * scalar surface. Each piece of work is a soft mound whose height is its mass
 * and whose width is fixed; where work clusters, the mounds sum and the terrain
 * rises. The result is a landscape whose shape is a measured property of the
 * corpus, not an invention.
 *
 * Pure. No canvas, no DOM, no clock. Given the same points it returns the same
 * field, which is what lets the contours be tested.
 */

/**
 * A Gaussian mound. Chosen over a linear or cubic falloff because the sum of
 * Gaussians is smooth everywhere — and contours extracted from a field with
 * discontinuous derivatives develop visible kinks exactly where mounds meet.
 *
 * @param {number} r2      Squared distance from the centre.
 * @param {number} sigma2  Squared width.
 */
function mound(r2, sigma2) {
  return Math.exp(-r2 / (2 * sigma2));
}

/**
 * Sample the field onto a regular grid.
 *
 * @param {Array<{x: number, y: number, mass: number}>} points  In [0,1]^2.
 * @param {object} opts
 * @param {number} opts.width       Grid columns.
 * @param {number} opts.height      Grid rows.
 * @param {number} [opts.sigma]     Mound width in field units. Larger merges
 *                                  clusters into landmasses; smaller isolates
 *                                  every point into its own island.
 * @param {Array<{x: number, y: number, weight: number}>} [opts.attractors]
 *                                  Extra positive mass — used to raise the
 *                                  terrain where the visitor's interest lies.
 * @returns {{data: Float32Array, width: number, height: number, min: number, max: number}}
 */
export function sampleField(points, opts) {
  const { width, height, sigma = 0.115, attractors = [] } = opts;
  if (!(width > 1 && height > 1)) throw new RangeError('sampleField: width and height must exceed 1');

  const sigma2 = sigma * sigma;
  const data = new Float32Array(width * height);
  let min = Infinity;
  let max = -Infinity;

  // Field units are normalised to the SHORTER axis so a wide canvas stretches
  // the terrain horizontally rather than distorting the mounds into ellipses.
  const aspect = width / height;

  for (let gy = 0; gy < height; gy += 1) {
    const fy = gy / (height - 1);
    for (let gx = 0; gx < width; gx += 1) {
      const fx = gx / (width - 1);
      let v = 0;

      for (let i = 0; i < points.length; i += 1) {
        const p = points[i];
        const dx = (fx - p.x) * aspect;
        const dy = fy - p.y;
        v += p.mass * mound(dx * dx + dy * dy, sigma2);
      }

      for (let i = 0; i < attractors.length; i += 1) {
        const a = attractors[i];
        const dx = (fx - a.x) * aspect;
        const dy = fy - a.y;
        // Attractors are wider and softer than work: they represent a region of
        // interest, not a single artefact.
        v += a.weight * mound(dx * dx + dy * dy, sigma2 * 2.6);
      }

      data[gy * width + gx] = v;
      if (v < min) min = v;
      if (v > max) max = v;
    }
  }

  return { data, width, height, min, max };
}

/**
 * Bilinear sample of a field at continuous coordinates.
 * @param {{data: Float32Array, width: number, height: number}} field
 */
export function sampleAt(field, x, y) {
  const { data, width, height } = field;
  const cx = Math.max(0, Math.min(width - 1.001, x));
  const cy = Math.max(0, Math.min(height - 1.001, y));
  const x0 = Math.floor(cx); const y0 = Math.floor(cy);
  const fx = cx - x0; const fy = cy - y0;
  const i = y0 * width + x0;
  const a = data[i]; const b = data[i + 1];
  const c = data[i + width]; const d = data[i + width + 1];
  return a * (1 - fx) * (1 - fy) + b * fx * (1 - fy) + c * (1 - fx) * fy + d * fx * fy;
}

/**
 * Gradient of the field, by central differences.
 * The contour direction is perpendicular to this, which is how strokes are made
 * to run ALONG the terrain rather than across it.
 * @returns {[number, number]}
 */
export function gradientAt(field, x, y) {
  const h = 1;
  return [
    (sampleAt(field, x + h, y) - sampleAt(field, x - h, y)) / (2 * h),
    (sampleAt(field, x, y + h) - sampleAt(field, x, y - h)) / (2 * h),
  ];
}