Description of the desired feature

Would be nice to implement the forward modelling method to compute the gravitational fields of tesseroids with variable densities described on Soler et al. (2019).

The main thing to do is adding the distance-based discretization algorithm, which should be implemented inside Numba jitted functions without allocating memory.
Also, the tesseroid_gravity function should take a density function per tesseroid as argument, which should be jitted in order to speed things up.

Are you willing to help implement and maintain this feature? Yes