/**
 * Interpolação bilinear 2D conforme plano técnico (sec. 3.2).
 *
 * Dados quatro pontos Q₁₁(x₁,y₁), Q₁₂(x₁,y₂), Q₂₁(x₂,y₁), Q₂₂(x₂,y₂),
 * estima f(x,y) por:
 *   f ≈ ((x₂-x)(y₂-y)·f₁₁ + (x-x₁)(y₂-y)·f₂₁ + (x₂-x)(y-y₁)·f₁₂ + (x-x₁)(y-y₁)·f₂₂) /
 *       ((x₂-x₁)(y₂-y₁))
 *
 * Aceita x fora do intervalo por extrapolação linear (clamp opcional).
 */

export type Grid2D = {
  xs: readonly number[];
  ys: readonly number[];
  values: readonly (readonly number[])[];
};

function findBracket(xs: readonly number[], x: number): [number, number, boolean] {
  const clamped = Math.max(xs[0], Math.min(x, xs[xs.length - 1]));
  const extrapolated = clamped !== x;
  if (xs.length === 1) return [0, 0, extrapolated];
  if (clamped >= xs[xs.length - 1]) {
    return [xs.length - 2, xs.length - 1, extrapolated];
  }
  for (let i = 0; i < xs.length - 1; i++) {
    const a = xs[i];
    const b = xs[i + 1];
    if (clamped >= a && clamped <= b) {
      return [i, i + 1, extrapolated];
    }
  }
  return [0, xs.length - 1, extrapolated];
}

export function bilinearInterp(grid: Grid2D, x: number, y: number): number {
  const { xs, ys, values } = grid;

  const [ix0, ix1] = findBracket(xs, x);
  const [iy0, iy1] = findBracket(ys, y);

  const x1 = xs[ix0];
  const x2 = xs[ix1];
  const y1 = ys[iy0];
  const y2 = ys[iy1];

  const f11 = values[iy0][ix0];
  const f21 = values[iy0][ix1];
  const f12 = values[iy1][ix0];
  const f22 = values[iy1][ix1];

  const dx = x2 - x1;
  const dy = y2 - y1;
  if (dx === 0 || dy === 0) return f11;

  const denom = dx * dy;
  const num =
    (x2 - x) * (y2 - y) * f11 +
    (x - x1) * (y2 - y) * f21 +
    (x2 - x) * (y - y1) * f12 +
    (x - x1) * (y - y1) * f22;

  return num / denom;
}