Efficient high-order accurate Fresnel diffraction via areal quadrature and the nonuniform FFT

We present a fast algorithm for computing the diffracted field from arbitrary binary (hard-edged) planar apertures and occulters in the scalar Fresnel approximation, for up to moderately high Fresnel numbers ($\lesssim 10^3$). It uses a high-order areal quadrature over the aperture, then exploits a single 2D nonuniform fast Fourier transform (NUFFT) to evaluate rapidly at target points (of order $10^7$ such points per second, independent of aperture complexity). It thus combines the high accuracy of edge integral methods with the high speed of Fourier methods. Its cost is ${\mathcal O}(n^2 \log n)$, where $n$ is the linear resolution required in source and target planes, to be compared with ${\mathcal O}(n^3)$ for edge integral methods. In tests with several aperture shapes, this translates to between 2 and 5 orders of magnitude acceleration. In starshade modeling for exoplanet astronomy, we find that it is roughly $10^4 \times$ faster than the state of the art in accurately computing the set of telescope pupil wavefronts. We provide a documented, tested MATLAB/Octave implementation. An appendix shows the mathematical equivalence of the boundary diffraction wave, angular integration, and line integral formulae, then analyzes a new non-singular reformulation that eliminates their common difficulties near the geometric shadow edge. This supplies a robust edge integral reference against which to validate the main proposal.