Smooth catmull quaternion interpolation

const EPSILON = 1e-6

// Rotation order: Yaw (Y) → Pitch (X) → Roll (Z)
function eulerToQuat(pitch, yaw, roll) {
    const d2r = Math.PI / 180

    // NEGATE yaw to match left-handed rotation around Y-axis
    const hy = -yaw * 0.5 * d2r
    const hp = pitch * 0.5 * d2r
    const hr = roll  * 0.5 * d2r

    // trig shortcuts
    const cy = Math.cos(hy), sy = Math.sin(hy)
    const cp = Math.cos(hp), sp = Math.sin(hp)
    const cr = Math.cos(hr), sr = Math.sin(hr)

    // q = qYaw * qPitch * qRoll
    const x = cy * sp * cr + sy * cp * sr
    const y = sy * cp * cr - cy * sp * sr
    const z = cy * cp * sr + sy * sp * cr
    const w = cy * cp * cr - sy * sp * sr

    return [x, y, z, w]
}

function normalizeQuat(q) {
    const mag = Math.sqrt(q.w * q.w + q.x * q.x + q.y * q.y + q.z * q.z)
    if (mag === 0) return { w: 1, x: 0, y: 0, z: 0 }
    return {
        w: q.w / mag,
        x: q.x / mag,
        y: q.y / mag,
        z: q.z / mag
    }
}

function slerpQuat(a, b, t) {
    let dot = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w

    if (dot < 0) {
        b = { x: -b.x, y: -b.y, z: -b.z, w: -b.w }
        dot = -dot
    }

    if (dot > 1 - EPSILON) {
        return normalizeQuat({
            x: a.x + (b.x - a.x) * t,
            y: a.y + (b.y - a.y) * t,
            z: a.z + (b.z - a.z) * t,
            w: a.w + (b.w - a.w) * t
        })
    }

    const theta_0 = Math.acos(dot)
    const sin_theta_0 = Math.sin(theta_0)

    const theta = theta_0 * t
    const sin_theta = Math.sin(theta)

    const s0 = Math.cos(theta) - dot * sin_theta / sin_theta_0
    const s1 = sin_theta / sin_theta_0

    return {
        x: s0 * a.x + s1 * b.x,
        y: s0 * a.y + s1 * b.y,
        z: s0 * a.z + s1 * b.z,
        w: s0 * a.w + s1 * b.w
    }
}

// Catmull-Rom spline interpolation for quaternions (w/x/y/z object format)
function catmull(q0, q1, q2, q3, t) {
    let a = slerpQuat(q0, q1, t + 1)
    let b = slerpQuat(q1, q2, t)
    let c = slerpQuat(q2, q3, t - 1)

    let d = slerpQuat(a, b, (t + 1) / 2)
    let e = slerpQuat(b, c, t / 2)

    return slerpQuat(d, e, t)
}