Add variable density tesseroids

[santisoler]

Add an implementation of the density-based discretization algorithm. Duplicate
the jit_tesseroid_gravity and the gauss_legendre_quadrature functions to
take a callable density argument. Reorganize the old and new tesseroid code,
moving some functions to private files. The dispatcher function chooses which
algorithm to use: constant density or variable density, serial or parallel. The
density-based discretization algorithm was modified to use the
scipy.optimize.minimize_scalar function to compute minimum and maximum
density values for normalization, and for computing the maximum absolute
difference. The density argument should be a function decorated with
numba.njit: we set this requirement to improve performance of the method. Add
example and test functions for the new feature. Add scipy as a dependency for
Harmonica.

Fixes #82

TODO:

• Add citation instructions for Soler2021 when using gradient-boosted equivalent sources

Reminders:

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

[santisoler]

I'm adding this comment in order to leave a record of experiences during this implementation.

Using the scipy.optimize.minimize_scalar function

I've experienced some errors when the minimize_scalar was trying to get the minimum and maximum value of a density within tesseroid's radial bounds, particularly with a strong L-shaped exponential function.
The function returned equal values for min and max, preventing the density-based discretization to split the tesseroid and impacting on the accuracy of the algorithm.
I solved it by comparing its outputs with the values at the boundaries and choose the minimum/maximum, depending on which value should be chosen.

Decorating the density function

The initial idea for this implementation was the tesseroid_gravity function to be able to get a density argument as a callable that will be njit decorated within the same function.
This behaviour has a strong impact on the performance of the algorithm: both the serial and parallel runs are considerably slower than the constant density case, and the parallel one is even slower than the serial one.

By asking a predecorated density function we solve this issue.

[santisoler]

@leouieda I think this is a good start. As always with tesseroids, I'm sure we will encounter some bugs or things to improve in the future. Let me know what do you think about the current status of this PR.

[leouieda]

@santisoler this looks really good! The tesseroid code is quite complex but there is no good way around it I guess. The only thing I noticed is the missing introductory text in the gallery example. That could be added in another PR as well to get this merged sooner.

[santisoler]

Thanks for the review @leouieda! I've added the introductory text and also changed the jit decorator for njit.

[leouieda]

Looks great! It's awesome to finally have this in here 🙂