// The cymbal design: a JavaScript script that returns an object describing
// the material, the bounding box, the outline, the thickness field, the
// stress-free profile, the mount, and the losses. The script runs in the page
// with a set of helper functions in scope (see HELPERS) and returns a plain
// object; compileDesign checks it and evaluatePlate samples the fields on a
// grid of the requested spacing.
//
// Lengths are meters throughout; `inch` and `mm` are available as constants.
// Every field function receives (x, y, r, theta): the point relative to the
// center of the bounding box, its distance from the center, and its angle
// from the x axis in radians.

import {
  type Grid,
  type InsideFn,
  type Loss,
  type Material,
  type Mount,
  type Plate,
  insideNode,
  makeGrid,
  nodeWeight,
  nodeX,
  nodeY,
} from "./params.ts";

export type PointFn = (x: number, y: number, r: number, theta: number) => number;
export type OutlineFn = (x: number, y: number, r: number, theta: number) => boolean;

/** A compiled design: the script's object, checked and with defaults filled in. */
export interface Design {
  name: string;
  material: Material;
  Lx: number;
  Ly: number;
  inside: OutlineFn | null;
  thickness: PointFn;
  profile: PointFn;
  mount: Mount | null;
  loss: Loss;
}

export const DEFAULT_LOSS: Loss = { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 };

// ---------- Helpers available to the script ----------

/** Deterministic uniform generator on [0, 1) (mulberry32). */
function mulberry32(seed: number): () => number {
  let a = Math.floor(seed) >>> 0;
  return () => {
    a = (a + 0x6d2b79f5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}

const NOISE_TERMS = 64;

/**
 * A smooth random field of unit RMS with the covariance exp(-d^2 / length^2):
 * a sum of plane waves whose wavevectors are Gaussian with standard deviation
 * sqrt(2) / length per component. The same seed gives the same field on any grid.
 */
export function noise2d(seed: number, length: number): (x: number, y: number) => number {
  const rnd = mulberry32(seed);
  const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());
  const L = Math.max(length, 1e-4);
  const terms = Array.from({ length: NOISE_TERMS }, () => ({ kx: (gauss() * Math.SQRT2) / L, ky: (gauss() * Math.SQRT2) / L, ph: 2 * Math.PI * rnd() }));
  const norm = 1 / Math.sqrt(NOISE_TERMS / 2);
  return (x, y) => {
    let v = 0;
    for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);
    return v * norm;
  };
}

export const HELPERS = {
  inch: 0.0254,
  mm: 1e-3,
  /** Outline of a circle of the given diameter, with an optional center hole of diameter `hole`. */
  circle: (diameter: number, hole = 0): OutlineFn => (_x, _y, r) => r <= diameter / 2 && r >= hole / 2,
  /** Outline of the ellipse with the given width and height, with an optional center hole. */
  ellipse: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => (x / (width / 2)) ** 2 + (y / (height / 2)) ** 2 <= 1 && r >= hole / 2,
  /** Outline of the rectangle of the given width and height (centered), with an optional center hole. */
  rectangle: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => Math.abs(x) <= width / 2 && Math.abs(y) <= height / 2 && r >= hole / 2,
  clamp: (x: number, lo: number, hi: number): number => Math.min(hi, Math.max(lo, x)),
  lerp: (a: number, b: number, t: number): number => a + (b - a) * t,
  /** 0 below e0, 1 above e1, a smooth cubic between. */
  smoothstep: (e0: number, e1: number, x: number): number => {
    const t = Math.min(1, Math.max(0, (x - e0) / (e1 - e0)));
    return t * t * (3 - 2 * t);
  },
  /** Gaussian bump of the given height, falling to 1/e at `radius` from the center; `d` is the distance from the center. */
  bump: (d: number, radius: number, height: number): number => height * Math.exp(-(d * d) / (radius * radius)),
  /** Cosine dome of the given height and base radius (zero slope at the rim and the apex); `d` is the distance from the center. */
  dome: (d: number, radius: number, height: number): number => (d >= radius ? 0 : 0.5 * height * (1 + Math.cos((Math.PI * d) / radius))),
  /** Parabolic bow of the given sag across a span L centered on 0; `s` is the coordinate along the span; the offset in -1..1 moves the crest toward one end. */
  bow: (s: number, L: number, sag: number, offset = 0): number => {
    const t = (2 * s) / L;
    return sag * (1 - t * t) * (1 + offset * t);
  },
  noise2d,
};

export type Helpers = typeof HELPERS;

// ---------- Compilation ----------

/** Message of an error thrown by the script or by the checks. */
export class DesignError extends Error {}

function num(v: unknown, what: string, min = -Infinity): number {
  if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} must be a finite number`);
  if (v < min) throw new DesignError(`${what} must be at least ${min}`);
  return v;
}

function fieldFn(v: unknown, what: string): PointFn {
  if (typeof v === "number") {
    if (!isFinite(v)) throw new DesignError(`${what} must be a finite number or a function`);
    return () => v;
  }
  if (typeof v === "function") return v as PointFn;
  throw new DesignError(`${what} must be a number or a function (x, y, r, theta) => number`);
}

/**
 * Runs the script and checks the object it returns. The script body has the
 * helpers in scope as free variables and must `return` the design object.
 */
export function compileDesign(source: string): Design {
  let raw: unknown;
  try {
    const fn = new Function(...Object.keys(HELPERS), source) as (...args: unknown[]) => unknown;
    raw = fn(...Object.values(HELPERS));
  } catch (e) {
    throw new DesignError(`script error: ${e instanceof Error ? e.message : String(e)}`);
  }
  if (!raw || typeof raw !== "object") throw new DesignError("the script must return an object (end it with `return { ... };`)");
  const o = raw as Record<string, unknown>;
  const name = typeof o.name === "string" && o.name.trim() ? o.name.trim() : "Untitled design";

  const m = o.material as Record<string, unknown> | undefined;
  if (!m || typeof m !== "object") throw new DesignError("material must be an object { E, nu, rho }");
  const material: Material = { E: num(m.E, "material.E", 1e6), nu: num(m.nu, "material.nu", 0), rho: num(m.rho, "material.rho", 1) };
  if (material.nu >= 0.5) throw new DesignError("material.nu must be below 0.5");

  const b = o.bounds as Record<string, unknown> | undefined;
  if (!b || typeof b !== "object") throw new DesignError("bounds must be an object { width, height }");
  const Lx = num(b.width, "bounds.width", 1e-3), Ly = num(b.height, "bounds.height", 1e-3);

  let inside: OutlineFn | null = null;
  if (o.inside !== undefined && o.inside !== null) {
    if (typeof o.inside !== "function") throw new DesignError("inside must be a function (x, y, r, theta) => boolean");
    inside = o.inside as OutlineFn;
  }

  if (o.thickness === undefined) throw new DesignError("thickness is required (a number in meters, or a function)");
  const thickness = fieldFn(o.thickness, "thickness");
  const profile = o.profile === undefined || o.profile === null ? () => 0 : fieldFn(o.profile, "profile");

  let mount: Mount | null = null;
  if (o.mount !== undefined && o.mount !== null && o.mount !== false) {
    const mo = o.mount as Record<string, unknown>;
    if (typeof mo !== "object") throw new DesignError("mount must be an object { x, y, radius, stiffness, damping } or null");
    mount = {
      x: mo.x === undefined ? 0 : num(mo.x, "mount.x"),
      y: mo.y === undefined ? 0 : num(mo.y, "mount.y"),
      radius: num(mo.radius, "mount.radius", 0),
      stiffness: num(mo.stiffness, "mount.stiffness", 0),
      damping: mo.damping === undefined ? 0 : num(mo.damping, "mount.damping", 0),
    };
  }

  const lo = (o.loss ?? {}) as Record<string, unknown>;
  if (typeof lo !== "object") throw new DesignError("loss must be an object { sigma0, sigma1, sigmaIn }");
  const loss: Loss = {
    sigma0: lo.sigma0 === undefined ? DEFAULT_LOSS.sigma0 : num(lo.sigma0, "loss.sigma0", 0),
    sigma1: lo.sigma1 === undefined ? DEFAULT_LOSS.sigma1 : num(lo.sigma1, "loss.sigma1", 0),
    sigmaIn: lo.sigmaIn === undefined ? DEFAULT_LOSS.sigmaIn : num(lo.sigmaIn, "loss.sigmaIn", 0),
  };

  return { name, material, Lx, Ly, inside, thickness, profile, mount, loss };
}

// ---------- Evaluation on a grid ----------

/** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
export function removePlane(f: Float64Array, g: Grid): void {
  let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
  for (let j = 0; j < g.ny; j++) {
    for (let i = 0; i < g.nx; i++) {
      const wgt = nodeWeight(g, i, j);
      if (wgt === 0) continue;
      const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
      sw += wgt; sx += wgt * x; sy += wgt * y;
      sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
      sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
    }
  }
  if (sw === 0) return;
  // Normal equations of the 3-parameter fit, solved by Cramer's rule.
  const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
  if (Math.abs(det) < 1e-300) return;
  const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
  const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
  const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) f[j * g.nx + i] -= a + b * nodeX(g, i) + c * nodeY(g, j);
}

function callField(fn: PointFn, x: number, y: number, what: string): number {
  let v: unknown;
  try {
    v = fn(x, y, Math.hypot(x, y), Math.atan2(y, x));
  } catch (e) {
    throw new DesignError(`${what} threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);
  }
  if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} returned ${String(v)} at (${x.toFixed(4)}, ${y.toFixed(4)})`);
  return v;
}

/**
 * Samples the design on a grid of about the requested spacing: the mask from
 * the outline, the thickness at the nodes (which must be positive on the
 * plate), and the profile at the nodes with its mean plane removed. Only the
 * slopes of the profile enter the physics, so its mean is immaterial, but its
 * tilt is not: the von Karman strains are not invariant under rotation, and a
 * tilted stress-free shape would add a spurious membrane coupling.
 */
export function evaluatePlate(d: Design, spacing: number): Plate {
  let inside: InsideFn | undefined;
  if (d.inside) {
    const fn = d.inside;
    inside = (x, y) => {
      try {
        return !!fn(x, y, Math.hypot(x, y), Math.atan2(y, x));
      } catch (e) {
        throw new DesignError(`inside threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);
      }
    };
  }
  const g = makeGrid(d.Lx, d.Ly, spacing, inside);
  const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
  let any = false;
  for (let j = 0; j < g.ny; j++) {
    for (let i = 0; i < g.nx; i++) {
      if (!insideNode(g, i, j)) continue;
      any = true;
      const x = nodeX(g, i), y = nodeY(g, j), k = j * g.nx + i;
      h[k] = callField(d.thickness, x, y, "thickness");
      if (h[k] <= 0) throw new DesignError(`thickness must be positive; it is ${h[k]} at (${x.toFixed(4)}, ${y.toFixed(4)})`);
      w0[k] = callField(d.profile, x, y, "profile");
    }
  }
  if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
  removePlane(w0, g);
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (!insideNode(g, i, j)) w0[j * g.nx + i] = 0;
  return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, mount: d.mount, loss: d.loss };
}

/** Peak-to-peak height and the largest slope of the profile over the plate nodes. */
export function profileStats(pl: Plate): { peakToPeak: number; maxSlope: number } {
  const g = pl.grid, w0 = pl.w0;
  let lo = Infinity, hi = -Infinity, maxSlope = 0;
  const at = (i: number, j: number) => w0[j * g.nx + i];
  const slope = (i: number, j: number, di: number, dj: number, hg: number) => {
    const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
    if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * hg);
    if (hasP) return (at(i + di, j + dj) - at(i, j)) / hg;
    if (hasM) return (at(i, j) - at(i - di, j - dj)) / hg;
    return 0;
  };
  for (let j = 0; j < g.ny; j++) {
    for (let i = 0; i < g.nx; i++) {
      if (!insideNode(g, i, j)) continue;
      const v = at(i, j);
      lo = Math.min(lo, v); hi = Math.max(hi, v);
      maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
    }
  }
  return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
}