// 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]];
}