The evaluation of highly oscillatory integrals is an important topic in many areas of computational wave propagation. The Numerical Steepest Descent (NSD) method is a powerful approach to computing such integrals, which combines complex contour deformation with quadrature rules. However, for fixed frequencies NSD may lose accuracy when stationary points are close with each other, or with endpoints of the integration contour. This issue can be dealt with standard quadrature inside neighbourhoods of stationary points, in which the number of oscillations is bounded and small, combined with NSD techniques away from stationary points. In this talk, I will present a simple “black-box” interface that automates contour deformation and integration, and describe a novel framework to rigorously analyse the numerical convergence of this class of methods.