Provided are techniques for numerically integrating an intensity distribution function over a numerical aperture in a manner dependent on a determination of whether the numerical aperture spans a Rayleigh singularity. Where a singularity exists, Gaussian quadrature (cubature) is performed using a set of weights and points (nodes) that account for the effect of the Wood anomaly present within the aperture space. The numerical aperture may be divided into subregions separated by curves where the Wood anomaly condition is satisfied. Each subregion is then numerically integrated and a weighted sum of the subregion contributions is the estimate of the integral. Alternatively, generalized Gaussian quadrature (cubature) is performed where an analytical polynomial function which accounts for the effect of the Wood anomaly present within the aperture space is integrated. Points and nodes generated from a fit of the analytical polynomial function are then used for integration of the intensity distribution function.