// Mirror of ../matlab/cymbal_profile.m: the cymbal meridian z(r), used by the
// 3D scene and the mode gallery. The physics reads the same geometry from the
// MATLAB side; this is only for drawing.
//
// Bow (r >= bell radius): a parabola in the rim-referred coordinate, flat at
// the rim. Bell (r < bell radius): the quartic that matches the bow's height,
// slope and curvature at the join and comes to a horizontal apex at r = 0.

export interface Meridian {
  rHole: number;
  rBell: number;
  rRim: number;
  z(r: number): number;
  thickness(r: number): number;
  apexHeight: number;
}

export function cymbalMeridian(
  radius: number,
  holeRadius: number,
  bellRadius: number,
  bellHeight: number,
  bowHeight: number,
  thicknessBell: number,
  thicknessRim: number,
): Meridian {
  const a = radius;
  const rh = Math.max(1e-4, Math.min(holeRadius, 0.5 * a));
  const rb = Math.min(Math.max(bellRadius, 2 * rh), 0.8 * a);
  const den = a - rb;
  const mb = (-2 * bowHeight) / den;
  const kb = (2 * bowHeight) / (den * den);
  // solve the 3x3 for the quartic coefficients (B r^2 + C r^3 + E r^4)
  const M = [
    [rb * rb, rb ** 3, rb ** 4],
    [2 * rb, 3 * rb * rb, 4 * rb ** 3],
    [2, 6 * rb, 12 * rb * rb],
  ];
  const rhs = [-bellHeight, mb, kb];
  const [B, C, E] = solve3(M, rhs);
  const A0 = bowHeight + bellHeight;
  return {
    rHole: rh,
    rBell: rb,
    rRim: a,
    apexHeight: A0,
    z(r: number): number {
      if (r < rb) return A0 + B * r * r + C * r ** 3 + E * r ** 4;
      const u = (a - r) / den;
      return bowHeight * u * u;
    },
    thickness(r: number): number {
      const taper = Math.min(1, Math.max(0, (a - r) / den));
      return thicknessRim + (thicknessBell - thicknessRim) * taper;
    },
  };
}

function solve3(M: number[][], b: number[]): [number, number, number] {
  const A = M.map((row, i) => [...row, b[i]]);
  for (let c = 0; c < 3; c++) {
    let piv = c;
    for (let r = c + 1; r < 3; r++) if (Math.abs(A[r][c]) > Math.abs(A[piv][c])) piv = r;
    [A[c], A[piv]] = [A[piv], A[c]];
    const d = A[c][c] || 1e-30;
    for (let k = c; k < 4; k++) A[c][k] /= d;
    for (let r = 0; r < 3; r++) {
      if (r === c) continue;
      const f = A[r][c];
      for (let k = c; k < 4; k++) A[r][k] -= f * A[c][k];
    }
  }
  return [A[0][3], A[1][3], A[2][3]];
}

export interface MeridianSamples {
  r: Float64Array;     // radius of each sample, m
  z: Float64Array;     // height above the rim plane, m
  sFrac: Float64Array; // arclength fraction from the hole, the model's strike_r
  total: number;       // meridian arclength, m
}

// Sample the meridian and tabulate arclength, so that the model's strike_r
// (an arclength fraction) can be converted to and from a radius. Both the 3D
// scene and the strike picker need this; the physics reads the same geometry
// from cymbal_profile.m.
export function meridianSamples(mer: Meridian, n = 256): MeridianSamples {
  const r = new Float64Array(n + 1);
  const z = new Float64Array(n + 1);
  const s = new Float64Array(n + 1);
  for (let i = 0; i <= n; i++) {
    r[i] = mer.rHole + (i / n) * (mer.rRim - mer.rHole);
    z[i] = mer.z(r[i]);
  }
  for (let i = 1; i <= n; i++) s[i] = s[i - 1] + Math.hypot(r[i] - r[i - 1], z[i] - z[i - 1]);
  const total = s[n];
  const sFrac = new Float64Array(n + 1);
  for (let i = 0; i <= n; i++) sFrac[i] = s[i] / total;
  return { r, z, sFrac, total };
}

// strike_r -> the point on the meridian, matching cymbal_point.m
export function pointAtFrac(m: MeridianSamples, frac: number) {
  const f = Math.min(0.999999, Math.max(0, frac));
  let i = 0;
  while (i < m.sFrac.length - 2 && m.sFrac[i + 1] < f) i++;
  const t = (f - m.sFrac[i]) / Math.max(1e-12, m.sFrac[i + 1] - m.sFrac[i]);
  const dr = m.r[i + 1] - m.r[i];
  const dz = m.z[i + 1] - m.z[i];
  const dl = Math.hypot(dr, dz) || 1;
  return { r: m.r[i] + t * dr, z: m.z[i] + t * dz, tr: dr / dl, tz: dz / dl };
}

// radius -> strike_r, the inverse the picker needs when the pointer moves
export function fracAtRadius(m: MeridianSamples, radius: number): number {
  const n = m.r.length - 1;
  const rr = Math.min(m.r[n], Math.max(m.r[0], radius));
  let i = 0;
  while (i < n - 1 && m.r[i + 1] < rr) i++;
  const t = (rr - m.r[i]) / Math.max(1e-12, m.r[i + 1] - m.r[i]);
  return m.sFrac[i] + t * (m.sFrac[i + 1] - m.sFrac[i]);
}