// Physical and numerical parameters of the simulation, plus derived quantities.

export interface PhysParams {
  Lx: number; // plate width along x [m]
  Ly: number; // plate height along y [m]
  h: number; // thickness [m]
  E: number; // Young's modulus [Pa]
  nu: number; // Poisson ratio
  rho: number; // density [kg/m^3]
  sigma0: number; // frequency-independent loss [1/s]
  sigma1: number; // frequency-dependent loss [m^2/s]
  sigmaIn: number; // additional loss of the in-plane motion [1/s]
  holeInset: number; // distance of the thread hole from the two edges of the top-left corner [m]
  pushAmount: number; // corner displacement [m]
  fingerRadius: number; // radius of the pushed patch [m]
  rampTime: number; // duration of the push ramp [s]
  holdMax: number; // maximum settle time before release [s]
  releaseTime: number; // time over which the finger force falls to zero [s]
  settleSigma0: number; // artificial damping during push and settle [1/s]
  settleKE: number; // kinetic energy threshold that counts as settled [J]
  spacing: number; // requested grid spacing [m]
  fs: number; // audio sample rate [Hz]
  fps: number; // movie frame rate (also the batch rate) [Hz]
  duration: number; // recorded duration after release [s]
  listenerDist: number; // [m]
  listenerAzimuth: number; // [deg], rotation about the plate's vertical axis (0 = on the normal)
  listenerElevation: number; // [deg]
  radiation: "rayleigh" | "dipole";
  bowShape: BowShape; // stress-free initial shape of the sheet
  bowAmplitude: number; // sag of the initial shape, center relative to the edges [m]
  excitation: Excitation;
  strikeFracX: number; // strike point as a fraction of the width from the left edge
  strikeFracY: number; // strike point as a fraction of the height from the bottom edge
  strikeVelocity: number; // stick tip speed at impact [m/s]
  stickMass: number; // effective mass of the stick [kg]
  contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
  contactExponent: number; // Hertz exponent alpha (1.5 for a sphere)
  contactRadius: number; // radius of the contact patch [m]
}

export type Excitation = "push" | "strike";

export type BowShape = "none" | "width" | "height" | "dish";

export const defaultParams: PhysParams = {
  Lx: 12 * 0.0254,
  Ly: 18 * 0.0254,
  h: 0.79e-3,
  E: 193e9,
  nu: 0.29,
  rho: 7900,
  sigma0: 1.15,
  sigma1: 5e-4,
  sigmaIn: 20,
  holeInset: 0.012,
  pushAmount: 0.01,
  fingerRadius: 0.01,
  rampTime: 0.1,
  holdMax: 0.5,
  releaseTime: 0.02,
  settleSigma0: 300,
  settleKE: 1e-8,
  spacing: 0.008,
  fs: 48000,
  fps: 120,
  duration: 4,
  listenerDist: 1.0,
  listenerAzimuth: 0,
  listenerElevation: 0,
  radiation: "rayleigh",
  bowShape: "width",
  bowAmplitude: 3e-3,
  excitation: "push",
  strikeFracX: 0.7,
  strikeFracY: 0.35,
  strikeVelocity: 3,
  stickMass: 0.03,
  contactStiffness: 5e8,
  contactExponent: 1.5,
  contactRadius: 5e-3,
};

/** Strike point in plate coordinates (relative to the center). */
export function strikePoint(p: PhysParams): [number, number] {
  return [(p.strikeFracX - 0.5) * p.Lx, (p.strikeFracY - 0.5) * p.Ly];
}

/** Stress-free initial shape w0 at every node. */
export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
  const w0 = new Float32Array(g.n);
  if (p.bowShape === "none" || p.bowAmplitude === 0) return w0;
  for (let j = 0; j < g.ny; j++) {
    for (let i = 0; i < g.nx; i++) {
      const X = (2 * nodeX(g, i)) / p.Lx; // -1..1
      const Y = (2 * nodeY(g, j)) / p.Ly;
      let f = 0;
      if (p.bowShape === "width") f = 1 - X * X;
      else if (p.bowShape === "height") f = 1 - Y * Y;
      else f = 1 - 0.5 * (X * X + Y * Y);
      w0[j * g.nx + i] = p.bowAmplitude * f;
    }
  }
  return w0;
}

export const SPEED_OF_SOUND = 343; // [m/s]
export const RHO_AIR = 1.2; // [kg/m^3]
export const HIST_LEN = 256; // ring buffer length for the retarded-time audio sum [samples]
export const RB_INTERVAL = 32; // steps between rigid-body projections while ringing
export const RB_INTERVAL_HELD = 4; // while the finger constraint is active

export interface Grid {
  nx: number;
  ny: number;
  hx: number;
  hy: number;
  n: number;
}

export interface Material {
  D: number; // bending stiffness [N m]
  C: number; // membrane stiffness E h / (1 - nu^2) [N/m]
  G: number; // membrane shear stiffness E h / (2 (1 + nu)) [N/m]
  nu: number;
  rhoh: number; // mass per unit area [kg/m^2]
}

export function makeGrid(p: PhysParams): Grid {
  const nx = Math.max(6, Math.round(p.Lx / p.spacing) + 1);
  const ny = Math.max(6, Math.round(p.Ly / p.spacing) + 1);
  return { nx, ny, hx: p.Lx / (nx - 1), hy: p.Ly / (ny - 1), n: nx * ny };
}

export function makeMaterial(p: PhysParams): Material {
  return {
    D: (p.E * p.h ** 3) / (12 * (1 - p.nu * p.nu)),
    C: (p.E * p.h) / (1 - p.nu * p.nu),
    G: (p.E * p.h) / (2 * (1 + p.nu)),
    nu: p.nu,
    rhoh: p.rho * p.h,
  };
}

/** Node coordinates relative to the plate center. */
export function nodeX(g: Grid, i: number): number {
  return (i - (g.nx - 1) / 2) * g.hx;
}
export function nodeY(g: Grid, j: number): number {
  return (j - (g.ny - 1) / 2) * g.hy;
}

/** Trapezoidal weight of a node: 1 in the interior, 1/2 on an edge, 1/4 at a corner. */
export function trapWeight(g: Grid, i: number, j: number): number {
  const tx = i === 0 || i === g.nx - 1 ? 0.5 : 1;
  const ty = j === 0 || j === g.ny - 1 ? 0.5 : 1;
  return tx * ty;
}

/** Decay time T60 [s] of a mode with angular frequency omega under the two-parameter loss model. */
export function t60(m: Material, p: PhysParams, omega: number): number {
  const kappa = Math.sqrt(m.D / m.rhoh);
  const k2 = omega / kappa;
  const sigma = p.sigma0 + p.sigma1 * k2;
  return (3 * Math.LN10) / sigma; // amplitude factor 1000 = 60 dB
}