A unified quadrature framework is developed for highly oscillatory integrals involving general \(n\) -fold products of Bessel functions on finite and semi-infinite intervals. The framework is intended for difficult regimes, including equal or nearly equal frequencies, algebraic endpoint behavior, and large Bessel orders. The method combines a left-end expansion, a selective Hankel splitting applied only to strongly oscillatory factors, and a unified steepest-descent quadrature operator. The operator uses Gauss–Legendre (GL) rules on short real segments and Gauss–Laguerre (GLag) rules on deformed rays. A split–merge strategy further reduces the number of subintegrals while preserving frequency awareness. The factor \(x^{\alpha}\) is handled directly, without desingularization. Separate bounds are given for the left-end truncation, finite-segment quadrature, vertical-ray quadrature, and tail contribution. These component bounds lead to an additive global error estimate and stability bounds for the composite scheme. Numerical experiments show exponential-type error decay with few nodes, robust behavior in equal-frequency and large-order regimes, and competitive cost–accuracy performance in the reported tests.
