// Numerical methods for transcendental equations

/**
 * Bisection method — finds root of f in [lo, hi] within tolerance.
 * Returns null if no sign change is found.
 */
export function bisect(f, lo, hi, tol = 1e-10, maxIter = 200) {
  let flo = f(lo)
  let fhi = f(hi)
  if (!isFinite(flo) || !isFinite(fhi)) return null
  if (flo * fhi > 0) return null // no sign change
  for (let i = 0; i < maxIter; i++) {
    const mid = (lo + hi) / 2
    if ((hi - lo) / 2 < tol) return mid
    const fmid = f(mid)
    if (fmid === 0) return mid
    if (flo * fmid < 0) { hi = mid; fhi = fmid }
    else { lo = mid; flo = fmid }
  }
  return (lo + hi) / 2
}

/**
 * Isentropic area-ratio function A/A* as a function of Mach M and gamma.
 *   A/A* = (1/M) * [(2/(γ+1)) * (1 + (γ-1)/2 * M²)]^((γ+1)/(2(γ-1)))
 */
export function areaRatioFromMach(M, gamma) {
  const exp = (gamma + 1) / (2 * (gamma - 1))
  const base = (2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M * M)
  return (1 / M) * Math.pow(base, exp)
}

/**
 * Solve Mach number from area ratio and gamma using bisection.
 * supersonic=true searches M > 1, false searches M < 1.
 * Returns null if unsolvable.
 */
export function machFromAreaRatio(eps, gamma, supersonic = true) {
  const f = (M) => areaRatioFromMach(M, gamma) - eps
  if (supersonic) return bisect(f, 1.0001, 200, 1e-9)
  else return bisect(f, 0.0001, 0.9999, 1e-9)
}