// Cymbal mid-surface height z0(x, y) and thickness h(x, y): bow, bell, rim warp and seeded
// hammer dents, as in cymsim/geometry.py. The slope that enters the strains is the CST
// gradient of the nodal heights (see shell.ts), so no analytic gradient is needed here.
import type { Imperfections, Profile } from "../config/schema.ts";
import { Rng } from "./rng.ts";

export class CymbalProfile {
  readonly dentX: Float64Array; readonly dentY: Float64Array; readonly dentAmp: Float64Array; readonly dentSigma: Float64Array;
  readonly warpM: Float64Array; readonly warpA: Float64Array; readonly warpPhase: Float64Array;
  readonly ovalPhase: number; readonly bellCx: number; readonly bellCy: number;

  readonly profile: Profile; readonly imperfections: Imperfections;

  constructor(profile: Profile, imperfections: Imperfections) {
    this.profile = profile; this.imperfections = imperfections;
    const p = profile, im = imperfections, rng = new Rng(im.seed);
    const r0 = p.bell_shoulder_radius + im.dent_margin, r1 = p.radius - im.dent_margin;
    const n = Math.max(0, Math.floor(im.n_dents));
    this.dentX = new Float64Array(n); this.dentY = new Float64Array(n);
    this.dentAmp = new Float64Array(n); this.dentSigma = new Float64Array(n);
    for (let i = 0; i < n; i++) {
      const rr = Math.sqrt(rng.uniform(r0 * r0, r1 * r1)), tt = rng.uniform(0, 2 * Math.PI);
      this.dentX[i] = rr * Math.cos(tt); this.dentY[i] = rr * Math.sin(tt);
    }
    for (let i = 0; i < n; i++) this.dentAmp[i] = im.dent_depth * rng.sign() * rng.uniform(0.5, 1.5);
    for (let i = 0; i < n; i++) this.dentSigma[i] = im.dent_sigma * rng.uniform(0.7, 1.3);
    const modes = im.rim_warp_modes.map((m) => Math.floor(m));
    this.warpM = Float64Array.from(modes);
    this.warpA = Float64Array.from(modes.map((m) => (im.rim_warp * 2) / m));
    this.warpPhase = Float64Array.from(modes.map(() => rng.uniform(0, 2 * Math.PI)));
    this.ovalPhase = rng.uniform(0, 2 * Math.PI);
    const bo = rng.uniform(0, 2 * Math.PI);
    this.bellCx = im.bell_offset * Math.cos(bo); this.bellCy = im.bell_offset * Math.sin(bo);
  }

  private shoulderRadius(phi: number): number {
    return this.profile.bell_shoulder_radius * (1 + this.imperfections.shoulder_ovality * Math.cos(2 * phi + this.ovalPhase));
  }

  /** 1 on the bell plateau, 0 beyond the shoulder, C1 cosine ramp in between. */
  private bellBlend(rBell: number, rSh: number): number {
    const top = this.profile.bell_top_radius;
    const t = Math.min(Math.max((rBell - top) / (rSh - top), 0), 1);
    return 0.5 * (1 + Math.cos(Math.PI * t));
  }

  private bowCoordinate(x: number, y: number) {
    const r = Math.sqrt(x * x + y * y + 1e-30);
    const phi = Math.atan2(y, r < 1e-6 ? 1e-6 : x);
    const rSh = this.shoulderRadius(phi);
    const rBell = Math.sqrt((x - this.bellCx) ** 2 + (y - this.bellCy) ** 2 + 1e-30);
    const s = Math.min(Math.max((this.profile.radius - r) / (this.profile.radius - rSh), 0), 1);
    return { s, phi, rBell, rSh };
  }

  /** Mid-surface height [m]. */
  z0(x: number, y: number): number {
    const p = this.profile;
    const { s, phi, rBell, rSh } = this.bowCoordinate(x, y);
    const u = Math.max(1 - s, 1e-12);
    let z = p.bow_height * (1 - Math.pow(u, p.bow_power)) + p.bell_height * this.bellBlend(rBell, rSh);
    let warp = 0;
    for (let k = 0; k < this.warpM.length; k++) warp += this.warpA[k] * Math.cos(this.warpM[k] * phi + this.warpPhase[k]);
    z += u * warp;
    for (let i = 0; i < this.dentX.length; i++) {
      const dx = x - this.dentX[i], dy = y - this.dentY[i], s2 = this.dentSigma[i] * this.dentSigma[i];
      z += this.dentAmp[i] * Math.exp(-(dx * dx + dy * dy) / (2 * s2));
    }
    return z;
  }

  /** Thickness [m]: linear taper along the bow, smooth step up into the bell. */
  thickness(x: number, y: number): number {
    const p = this.profile;
    const { s, rBell, rSh } = this.bowCoordinate(x, y);
    return p.h_rim + (p.h_shoulder - p.h_rim) * s + (p.h_bell - p.h_shoulder) * this.bellBlend(rBell, rSh);
  }
}