Improved Fourier modeling of gravity fields caused by polyhedral bodies: with applications to asteroid Bennu and comet 67P/Churyumov–Gerasimenko

This paper presents an improved algorithm for the 2D and 3D Fourier forward modeling of gravity fields caused by polyhedral bodies with constant and exponential density distributions. Three modifications have been made to the Fourier forward algorithm introduced in a previous paper. First, vertex-based Fourier-domain expressions are used instead of the original face-based Fourier-domain expressions, which simplify the computation of the anomaly spectrum considerably, especially in 3D modeling problems. Second, instead of using a pure Gauss-FFT sampling of the anomaly spectrum, we apply an improved sampling strategy by combining a nonuniform spherical sampling with a low-order Gauss-FFT sampling. In this way, the number of samplings required in the Fourier domain reduces to about \(\frac{1}{3}\) and \(\frac{1}{7}\) of those required in a pure Gauss-FFT algorithm for 2D and 3D modeling problems, respectively. A significant acceleration over the original algorithm is achieved. Third, we incorporate all three types of nonuniform fast Fourier transform algorithms to transform directly a uniform or nonuniform anomaly spectrum to gravity fields either on a regular grid, or at a set of arbitrary positions. Extra interpolation operations are no longer needed. Synthetic numerical tests show that for gravity vector components, the new algorithm runs about 3 times faster in 2D modeling and 7 times faster in 3D modeling than the original ones, while maintaining the same level of accuracy. For the gravity potential, the new algorithm is significantly superior to the pure Gauss-FFT solution both in numerical accuracy and in efficiency. We apply this novel approach to compute the gravitational fields of asteroid 101955 Bennu and comet \(67\hbox {P/Churyumov}\)–Gerasimenko. The 2D algorithm works very efficiently for the computation of gravity fields on horizontal planes. The 3D algorithm is valid both outside, on, and inside the source’s bounding surface, with relative errors less than 0.1% for the gravity potential and less than 2% for the gravity vector. By comparing to modeling results of analytical and spherical harmonic-based solutions, we generally conclude that the Fourier-based algorithm introduced here is an attractive alternative to these conventional solutions, especially for nonspherical, irregularly shaped bodies with complex geometries.