// Stress-free initial shape of the sheet. The shape is a list of features
// (bows, a twist, dents, creases, a random field) summed at the nodes, with
// the mean plane removed. Only the slopes of the shape enter the physics, so
// the mean is immaterial; the tilt is not, because the von Karman strains are
// not invariant under rotation, and a tilted stress-free shape would add a
// spurious membrane coupling. Every shape is therefore made plane-free here.

import { type Grid, type PhysParams, insideNode, nodeWeight, nodeX, nodeY } from "./params.ts";

export type FeatureType = "bow" | "twist" | "dent" | "crease" | "noise";

/** Parabolic bow across the width (axis x) or the height (axis y); the offset in -1..1 moves the crest toward one edge. */
export interface BowFeature { type: "bow"; axis: "x" | "y"; sag: number; offset: number }
/** Saddle x*y with the given height at the corners (alternating in sign). */
export interface TwistFeature { type: "twist"; height: number }
/** Gaussian bump at (x, y) given as fractions of the width and the height, with the 1/e radius and the height. */
export interface DentFeature { type: "dent"; x: number; y: number; radius: number; height: number }
/** Fold along the line through (x, y) at `angle` degrees from the x axis, kinked by `fold` degrees (positive: ridge), rounded over `width`. */
export interface CreaseFeature { type: "crease"; x: number; y: number; angle: number; fold: number; width: number }
/** Seeded random field with Gaussian correlation of length `length` and the given RMS over the sheet. */
export interface NoiseFeature { type: "noise"; seed: number; rms: number; length: number }
export type ShapeFeature = BowFeature | TwistFeature | DentFeature | CreaseFeature | NoiseFeature;

export interface FeatureField {
  key: string;
  label: string;
  unit: string;
  scale: number; // form value = SI value / scale
  step: number;
  options?: [string, string][]; // a choice rather than a number
}

export interface FeatureSpec {
  label: string;
  description: string;
  fields: FeatureField[];
}

/** Presentation and encoding order of the fields of each feature type. */
export const FEATURE_SPECS: Record<FeatureType, FeatureSpec> = {
  bow: {
    label: "Bow",
    description: "Parabolic bow with zero sag along two opposite edges. The apex offset (-1 to 1) moves the crest toward one edge.",
    fields: [
      { key: "axis", label: "Across", unit: "", scale: 1, step: 1, options: [["x", "the width"], ["y", "the height"]] },
      { key: "sag", label: "Sag", unit: "mm", scale: 1e-3, step: 0.5 },
      { key: "offset", label: "Apex offset", unit: "", scale: 1, step: 0.1 },
    ],
  },
  twist: {
    label: "Twist",
    description: "Saddle shape: two opposite corners up by the height, the other two down.",
    fields: [{ key: "height", label: "Height at the corners", unit: "mm", scale: 1e-3, step: 0.5 }],
  },
  dent: {
    label: "Dent",
    description: "Gaussian bump; the radius is where the height has fallen to 1/e. A negative height is a dent toward the back.",
    fields: [
      { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
      { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
      { key: "radius", label: "Radius", unit: "mm", scale: 1e-3, step: 5 },
      { key: "height", label: "Height", unit: "mm", scale: 1e-3, step: 0.5 },
    ],
  },
  crease: {
    label: "Crease",
    description: "Fold along a line: the two sides tilt by half the fold angle each. A positive fold is a ridge, a negative one a valley; the width rounds the kink.",
    fields: [
      { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
      { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
      { key: "angle", label: "Line angle from x axis", unit: "deg", scale: 1, step: 15 },
      { key: "fold", label: "Fold angle", unit: "deg", scale: 1, step: 0.5 },
      { key: "width", label: "Rounding width", unit: "mm", scale: 1e-3, step: 5 },
    ],
  },
  noise: {
    label: "Random field",
    description: "Smooth random undulation with Gaussian correlation; the same seed gives the same field on any grid.",
    fields: [
      { key: "seed", label: "Seed", unit: "", scale: 1, step: 1 },
      { key: "rms", label: "RMS amplitude", unit: "mm", scale: 1e-3, step: 0.1 },
      { key: "length", label: "Correlation length", unit: "mm", scale: 1e-3, step: 10 },
    ],
  },
};

export const FEATURE_TYPES = Object.keys(FEATURE_SPECS) as FeatureType[];

/** A feature of the given type with sensible default values. */
export function newFeature(type: FeatureType): ShapeFeature {
  switch (type) {
    case "bow": return { type, axis: "x", sag: 3e-3, offset: 0 };
    case "twist": return { type, height: 1e-3 };
    case "dent": return { type, x: 0.5, y: 0.5, radius: 0.03, height: 1e-3 };
    case "crease": return { type, x: 0.5, y: 0.5, angle: 30, fold: 1, width: 0.015 };
    case "noise": return { type, seed: 1, rms: 0.5e-3, length: 0.06 };
  }
}

export function cloneFeatures(features: ShapeFeature[]): ShapeFeature[] {
  return features.map((f) => ({ ...f }));
}

/** The shapes of the earlier interface: a parabolic bow across the width, the height, or both (dish). */
export function legacyFeatures(bowShape: string, amplitude: number): ShapeFeature[] {
  if (bowShape === "none" || amplitude === 0) return [];
  if (bowShape === "height") return [{ type: "bow", axis: "y", sag: amplitude, offset: 0 }];
  if (bowShape === "dish") return [{ type: "bow", axis: "x", sag: amplitude / 2, offset: 0 }, { type: "bow", axis: "y", sag: amplitude / 2, offset: 0 }];
  return [{ type: "bow", axis: "x", sag: amplitude, offset: 0 }];
}

// ---------- Text encoding (URL) ----------

function fmt(v: number): string {
  return v === 0 ? "0" : String(+v.toPrecision(8));
}

/** Compact text form: `type:v1,v2,...;type:...` with the values in SI, or `none` for a flat sheet. */
export function encodeFeatures(features: ShapeFeature[]): string {
  if (features.length === 0) return "none";
  return features
    .map((f) => `${f.type}:${FEATURE_SPECS[f.type].fields.map((fld) => {
      const v = (f as unknown as Record<string, number | string>)[fld.key];
      return typeof v === "number" ? fmt(v) : String(v);
    }).join(",")}`)
    .join(";");
}

/** Inverse of encodeFeatures. Returns null when nothing in the text can be read; entries that cannot be read are skipped. */
export function parseFeatures(text: string): ShapeFeature[] | null {
  const s = text.trim();
  if (s === "" || s === "none" || s === "flat") return [];
  const out: ShapeFeature[] = [];
  let any = false;
  for (const item of s.split(";")) {
    const colon = item.indexOf(":");
    if (colon < 0) continue;
    const type = item.slice(0, colon).trim() as FeatureType;
    const spec = FEATURE_SPECS[type];
    if (!spec) continue;
    const vals = item.slice(colon + 1).split(",").map((v) => v.trim());
    if (vals.length !== spec.fields.length) continue;
    const f = newFeature(type) as unknown as Record<string, number | string>;
    let ok = true;
    spec.fields.forEach((fld, i) => {
      if (fld.options) {
        if (!fld.options.some(([v]) => v === vals[i])) ok = false;
        else f[fld.key] = vals[i];
      } else {
        const num = parseFloat(vals[i]);
        if (!isFinite(num)) ok = false;
        else f[fld.key] = num;
      }
    });
    if (!ok) continue;
    out.push(f as unknown as ShapeFeature);
    any = true;
  }
  return any ? out : null;
}

// ---------- Evaluation ----------

/** 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;

/**
 * Adds the field of one feature to f (float64 per node, in meters). Nodes
 * outside the sheet are filled too and zeroed at the end.
 */
function addFeature(f: Float64Array, feat: ShapeFeature, p: PhysParams, g: Grid): void {
  const { nx, ny } = g;
  if (feat.type === "noise") {
    // Sum of random plane waves whose wavevectors are Gaussian with standard
    // deviation sqrt(2)/L per component, which gives the covariance exp(-r^2/L^2).
    const rnd = mulberry32(feat.seed);
    const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());
    const L = Math.max(feat.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 raw = new Float64Array(g.n);
    let sw = 0, s1 = 0, s2 = 0;
    for (let j = 0; j < ny; j++) {
      for (let i = 0; i < nx; i++) {
        const x = nodeX(g, i), y = nodeY(g, j);
        let v = 0;
        for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);
        raw[j * nx + i] = v;
        const wgt = nodeWeight(g, i, j);
        sw += wgt; s1 += wgt * v; s2 += wgt * v * v;
      }
    }
    const mean = s1 / sw, rms = Math.sqrt(Math.max(0, s2 / sw - mean * mean));
    const scale = rms > 0 ? feat.rms / rms : 0;
    for (let k = 0; k < g.n; k++) f[k] += scale * (raw[k] - mean);
    return;
  }
  for (let j = 0; j < ny; j++) {
    for (let i = 0; i < nx; i++) {
      const x = nodeX(g, i), y = nodeY(g, j);
      const X = (2 * x) / p.Lx, Y = (2 * y) / p.Ly; // -1..1 over the bounding box
      let v = 0;
      switch (feat.type) {
        case "bow": {
          const s = feat.axis === "x" ? X : Y;
          v = feat.sag * (1 - s * s) * (1 + feat.offset * s);
          break;
        }
        case "twist":
          v = feat.height * X * Y;
          break;
        case "dent": {
          const dx = x - (feat.x - 0.5) * p.Lx, dy = y - (feat.y - 0.5) * p.Ly;
          const r = Math.max(feat.radius, 1e-4);
          v = feat.height * Math.exp(-(dx * dx + dy * dy) / (r * r));
          break;
        }
        case "crease": {
          const th = (feat.angle * Math.PI) / 180;
          const d = -(x - (feat.x - 0.5) * p.Lx) * Math.sin(th) + (y - (feat.y - 0.5) * p.Ly) * Math.cos(th);
          const s = Math.max(feat.width, 0);
          v = -Math.tan((feat.fold * Math.PI) / 360) * (Math.sqrt(d * d + s * s) - s);
          break;
        }
      }
      f[j * nx + i] += v;
    }
  }
}

/** Subtracts the mass-weighted least-squares plane a + b x + c y from f over the sheet 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);
}

/** Stress-free initial shape w0 at every node: the sum of the features with the mean plane removed; zero outside the sheet. */
export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
  const w0 = new Float32Array(g.n);
  if (p.features.length === 0) return w0;
  const f = new Float64Array(g.n);
  for (const feat of p.features) addFeature(f, feat, p, g);
  removePlane(f, 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] = f[j * g.nx + i];
  return w0;
}

/** Peak-to-peak height and the largest slope of a shape over the sheet nodes. */
export function shapeStats(w0: ArrayLike<number>, g: Grid): { peakToPeak: number; maxSlope: number } {
  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, h: 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 * h);
    if (hasP) return (at(i + di, j + dj) - at(i, j)) / h;
    if (hasM) return (at(i, j) - at(i - di, j - dj)) / h;
    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 };
}

// ---------- Random configurations ----------

/**
 * A random feature list whose heights are of the order of `scale` (m): always
 * a bow across one axis, sometimes one across the other, a twist, up to two
 * dents, a crease, and a random field. `rnd` is a uniform generator on [0, 1).
 */
export function randomFeatures(p: PhysParams, scale: number, rnd: () => number = Math.random): ShapeFeature[] {
  const U = (a: number, b: number) => a + (b - a) * rnd();
  const sign = () => (rnd() < 0.5 ? -1 : 1);
  const Lmin = Math.min(p.Lx, p.Ly);
  const out: ShapeFeature[] = [];
  const first: "x" | "y" = rnd() < 0.5 ? "x" : "y";
  out.push({ type: "bow", axis: first, sag: sign() * U(0.3, 1) * scale, offset: U(-0.6, 0.6) });
  if (rnd() < 0.5) out.push({ type: "bow", axis: first === "x" ? "y" : "x", sag: sign() * U(0.1, 0.5) * scale, offset: U(-0.6, 0.6) });
  if (rnd() < 0.5) out.push({ type: "twist", height: sign() * U(0.1, 0.5) * scale });
  const dents = rnd() < 0.3 ? 0 : rnd() < 0.6 ? 1 : 2;
  for (let k = 0; k < dents; k++) out.push({ type: "dent", x: U(0.15, 0.85), y: U(0.15, 0.85), radius: U(0.1, 0.2) * Lmin, height: sign() * U(0.2, 0.7) * scale });
  if (rnd() < 0.3) {
    // A rise of a fraction of the scale over a quarter of the smaller side.
    const slope = (U(0.2, 0.6) * scale) / (Lmin / 4);
    out.push({ type: "crease", x: U(0.3, 0.7), y: U(0.3, 0.7), angle: U(0, 180), fold: sign() * ((2 * Math.atan(slope) * 180) / Math.PI), width: 0.05 * Lmin });
  }
  if (rnd() < 0.7) out.push({ type: "noise", seed: 1 + Math.floor(rnd() * 99999), rms: U(0.05, 0.25) * scale, length: U(0.1, 0.35) * Lmin });
  // Round the values to what the form shows, so the list reads the same in the form and in a link.
  for (const f of out) {
    const rec = f as unknown as Record<string, number | string>;
    for (const fld of FEATURE_SPECS[f.type].fields) {
      if (fld.options || fld.key === "seed") continue;
      const v = rec[fld.key] as number;
      rec[fld.key] = +((v / fld.scale).toPrecision(4)) * fld.scale;
    }
  }
  return out;
}