// In-browser check of the WGSL step kernel against the float64 reference:
// one step from a smooth random state of a few thicknesses amplitude, with
// the losses and the stand (lift-off included), compared as accelerations.

import { type Plate, type SimParams, meanThickness, nodeX, nodeY, insideNode, standRates } from "../params.ts";
import { accelerations, addDamping, fields, fullLoss } from "../model/reference.ts";
import { fitSigma, decayRate } from "../loss.ts";
import { Simulation } from "./simulation.ts";

export interface SelfTestResult {
  n: number;
  /** Largest acceleration error over all nodes and components relative to the largest reference acceleration. */
  maxRel: number;
  dt: number;
}

/** A smooth random field of the given amplitude on the plate nodes (zero outside). */
export function smoothField(pl: Plate, amplitude: number, seed: number): Float64Array {
  const g = pl.grid;
  const f = new Float64Array(g.n);
  let s = seed >>> 0;
  const rnd = () => { s = (s * 1664525 + 1013904223) >>> 0; return s / 4294967296; };
  const terms = Array.from({ length: 6 }, () => ({ kx: (rnd() - 0.5) * 60, ky: (rnd() - 0.5) * 60, ph: rnd() * 2 * Math.PI, a: rnd() }));
  for (let j = 0; j < g.ny; j++) {
    for (let i = 0; i < g.nx; i++) {
      if (!insideNode(g, i, j)) continue;
      const x = nodeX(g, i), y = nodeY(g, j);
      let v = 0;
      for (const t of terms) v += t.a * Math.sin(t.kx * x + t.ky * y + t.ph);
      f[j * g.nx + i] = v * amplitude / 3;
    }
  }
  return f;
}

export async function runSelfTest(device: GPUDevice, pl: Plate, params: SimParams): Promise<SelfTestResult> {
  const g = pl.grid, n = g.n;
  const hm = meanThickness(pl);
  const fit = fitSigma(pl.loss, pl.mat, hm);
  const sr = standRates(pl);
  const loss = fullLoss(fit, sr, decayRate(pl.loss, 5000));
  const sim = new Simulation(device, pl, params, loss);
  const amp = 2.5 * hm;
  const f = { u: smoothField(pl, 0.02 * amp, 1), v: smoothField(pl, 0.02 * amp, 2), w: smoothField(pl, amp, 3) };
  const vel = { u: smoothField(pl, 0.05, 4), v: smoothField(pl, 0.05, 5), w: smoothField(pl, 1.0, 6) };
  const pos32 = new Float32Array(4 * n), vel32 = new Float32Array(3 * n);
  for (let k = 0; k < n; k++) {
    pos32[k] = f.u[k]; pos32[n + k] = f.v[k]; pos32[2 * n + k] = f.w[k];
    vel32[k] = vel.u[k]; vel32[n + k] = vel.v[k]; vel32[2 * n + k] = vel.w[k];
  }
  // Use the float32-rounded state as the reference input.
  for (let k = 0; k < n; k++) {
    f.u[k] = pos32[k]; f.v[k] = pos32[n + k]; f.w[k] = pos32[2 * n + k];
    vel.u[k] = vel32[k]; vel.v[k] = vel32[n + k]; vel.w[k] = vel32[2 * n + k];
  }
  sim.writeState(pos32, vel32);
  await sim.stepOnce();
  const st = await sim.readState();
  const ref = fields(n);
  accelerations(pl, f, ref, params.padLiftoff, sr);
  addDamping(pl, loss, vel, ref, sr);
  let maxRef = 0, maxErr = 0;
  const dt = sim.dt;
  for (let k = 0; k < n; k++) {
    const au = (st.vel[k] - vel.u[k]) / dt, av = (st.vel[n + k] - vel.v[k]) / dt, aw = (st.vel[2 * n + k] - vel.w[k]) / dt;
    maxRef = Math.max(maxRef, Math.abs(ref.u[k]), Math.abs(ref.v[k]), Math.abs(ref.w[k]));
    maxErr = Math.max(maxErr, Math.abs(au - ref.u[k]), Math.abs(av - ref.v[k]), Math.abs(aw - ref.w[k]));
  }
  sim.destroy();
  return { n, maxRel: maxErr / maxRef, dt };
}