Based on a brief review of forward algorithms for the computation of topographic gravitational and magnetic effects, including spatial, spectral and hybrid-domain algorithms working in either Cartesian or spherical coordinate systems, we introduce a new algorithm, namely the CP-FFT algorithm, for fast computation of terrain-induced gravitational and magnetic effects on arbitrary undulating surfaces. The CP-FFT algorithm, working in the hybrid spatial-spectral domain, is based on a combination of CANDECOMP/PARAFAC (CP) tensor decomposition of gravitational integral kernels and 2D Fast Fourier Transform (FFT) evaluation of discrete convolutions. By replacing the binomial expansion in classical FFT-based terrain correction algorithms using CP decomposition, convergence of the outer-zone computation can be achieved with significantly reduced inner-zone radius. Additionally, a Gaussian quadrature mass line model is introduced to accelerate the computation of the inner zone effect. We validate our algorithm by computing the gravitational potential, the gravitational vector, the gravity gradient tensor, and magnetic fields caused by densely-sampled topographic and bathymetric digital elevation models of selected mountainous areas around the globe. Both constant and variable density/magnetization models, with computation surfaces on, above and below the topography are considered. Comparisons between our new method and space-domain rigorous solutions show that with modeling errors well below existing instrumentation error levels, the calculation speed is accelerated thousands of times in all numerical tests. We release a set of open-source code written in MATLAB language to meet the needs of geodesists and geophysicists in related fields to carry out more efficiently topographic modeling in Cartesian coordinates under planar approximation.
