  Paraxial diffraction modeling based on the Fourier transform has seen
widespread implementation for simulating the response of a diffraction-limited
optical system. For systems where the paraxial assumption is not sufficient, a
class of algorithms has been developed that employs hybrid propagation physics
to compute the propagation of an elementary beamlet along geometric ray paths.
These "beamlet decomposition" algorithms include the well-known Gaussian
Beamlet Decomposition (GBD) algorithm, of which several variants have been
created. To increase the computational efficiency of the GBD algorithm, we
derive an alternative expression of the technique that utilizes the analytical
propagation of beamlets to tilted planes. We then use this accelerated
algorithm to conduct a parameter-space search to find the optimal combination
of free parameters in GBD to construct the analytical Airy function. The
experiment is conducted on a consumer-grade CPU, and a high-performance GPU,
where the new algorithm is 34 times faster than the previously published
algorithm on CPUs, and 67,513 times faster on GPUs.
