Forward modelling for point masses in Cartesian coordiantes (#71)

Extend existing point mass forward modelling to Cartesian coordinates. Improve
docstring of point_mass_gravity function by adding mathematical definitions and
references. Add tests and a gallery example for the extended feature.

doc/references.rst

@@ -3,7 +3,9 @@ References
 
 .. [AmanteEakins2009] Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:`10.7289/V5C8276M <https://doi.org/10.7289/V5C8276M>`__
 .. [BarthelmesKohler2016] Barthelmes, F. and Kohler, W. (2016), International Centre for Global Earth Models (ICGEM), in: Drewes, H., Kuglitsch, F., Adam, J. et al., The Geodesists Handbook 2016, Journal of Geodesy (2016), 90(10), pp 907-1205, doi:`10.1007/s00190-016-0948-z <https://doi.org/10.1007/s00190-016-0948-z>`__
+.. [Blakely1995] Blakely, R. (1995). Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. doi:`10.1017/CBO9780511549816 <https://doi.org/10.1017/CBO9780511549816>`__
 .. [Forste_etal2014] Förste, Christoph; Bruinsma, Sean.L.; Abrikosov, Oleg; Lemoine, Jean-Michel; Marty, Jean Charles; Flechtner, Frank; Balmino, G.; Barthelmes, F.; Biancale, R. (2014): EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse. GFZ Data Services. doi:`10.5880/icgem.2015.1 <http://doi.org/10.5880/icgem.2015.1>`__
+.. [Grombein2013] Grombein, T., Seitz, K., Heck, B. (2013), Optimized formulas for the gravitational field of a tesseroid, Journal of Geodesy. doi:`10.1007/s00190-013-0636-1 <https://doi.org/10.1007/s00190-013-0636-1>`__
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [LiGotze2001] Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi:`10.1190/1.1487109 <https://doi.org/10.1190/1.1487109>`__
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.

examples/forward/point_mass.py

@@ -0,0 +1,54 @@
+"""
+Point Mass in Cartesian Coordinates
+===================================
+
+The simplest geometry used to compute gravitational fields are point masses.
+Although modelling geologic structures with point masses can be challenging, they are
+very usefull for other purposes: creating synthetic models, solving inverse problems
+very quickly, generating equivalent sources for interpolation or gridding, etc.
+The gravitational fields generated by point masses can be quickly computed
+either in Cartesian or in geocentric spherical coordinate systems.
+We will compute the downward component of the gravitational acceleration generated by
+a set of point masses on a computation grid given in Cartesian coordinates. We will do
+it throught the :func:`harmonica.point_mass_gravity` function.
+"""
+import harmonica as hm
+import verde as vd
+import matplotlib.pyplot as plt
+import matplotlib.ticker
+
+
+# Define two point masses with oposite mass values of 10000000 kg
+easting = [5e3, 15e3]
+northing = [5e3, 15e3]
+down = [7e3, 2.5e3]
+points = [easting, northing, down]
+# We're using "negative" masses to represent a "mass deficit"
+masses = [10e6, -10e6]
+
+# Define computation points on a grid at 500m above the ground
+coordinates = vd.grid_coordinates(
+ region=[0, 20e3, 0, 20e3], shape=(80, 80), extra_coords=-500
+)
+
+# Compute the downward component of the acceleration
+gravity = hm.point_mass_gravity(
+ coordinates, points, masses, field="g_z", coordinate_system="cartesian"
+)
+print(gravity)
+
+# Plot the gravitational field
+fig, ax = plt.subplots(figsize=(8, 9))
+ax.set_aspect("equal")
+img = plt.pcolormesh(*coordinates[:2], gravity)
+# Add colorbar
+fmt = matplotlib.ticker.ScalarFormatter(useMathText=True)
+fmt.set_powerlimits((0, 0))
+plt.colorbar(img, ax=ax, format=fmt, pad=0.04, shrink=0.73, label="mGal")
+ax.set_title("Downward component of gravitational acceleration")
+# Convert axes units to km
+ax.set_xticklabels(ax.get_xticks() * 1e-3)
+ax.set_yticklabels(ax.get_yticks() * 1e-3)
+ax.set_xlabel("km")
+ax.set_ylabel("km")
+plt.show()

harmonica/forward/point_mass.py

@@ -7,28 +7,117 @@
 from ..constants import GRAVITATIONAL_CONST
 
 
-def point_mass_gravity(coordinates, points, masses, field, dtype="float64"):
- """
- Compute gravitational fields of a point mass on spherical coordinates.
+def point_mass_gravity(
+ coordinates, points, masses, field, coordinate_system="cartesian", dtype="float64"
+):
+ r"""
+ Compute gravitational fields of point masses.
+
+ It can compute the gravitational fields of point masses on a set of computation
+ points defined either in Cartesian or geocentric spherical coordinates.
+
+ The potential gravity field generated by a point mass with mass :math:`m` located at
+ a point :math:`Q` on a computation point :math:`P` can be computed as:
+
+ .. math::
+
+ V(P) = \frac{G m}{l},
+
+ where :math:`G` is the gravitational constant and :math:`l` is the Euclidean
+ distance between :math:`P` and :math:`Q` [Blakely1995]_.
+
+ In Cartesian coordinates, the points :math:`P` and :math:`Q` are given by :math:`x`,
+ :math:`y` and :math:`z` coordinates, which can be translated into ``northing``,
+ ``easting`` and ``down``, respectively.
+ If :math:`P` is located at :math:`(x, y, z)`, and :math:`Q` at :math:`(x_p, y_p,
+ z_p)`, the distance :math:`l` can be computed as:
+
+ .. math::
+
+ l = \sqrt{ (x - x_p)^2 + (y - y_p)^2 + (z - z_p)^2 }.
+
+ The gradient of the potential, also known as the gravity acceleration vector
+ :math:`\vec{g}`, is defined as:
+
+ .. math::
+
+ \vec{g} = \nabla V.
+
+ Therefore, the :math:`z` component of :math:`\vec{g}` at the point :math:`P` can be
+ computed as (remember that :math:`z` points downward):
+
+ .. math::
+
+ g_z(P) = \frac{G m}{l^3} (z_p - z).
+
+ On a geocentric spherical coordinate system, the points :math:`P` and :math:`Q` are
+ given by the ``longitude``, ``latitude`` and ``radius`` coordinates, i.e.
+ :math:`\lambda`, :math:`\varphi` and :math:`r`, respectively. On this coordinate
+ system, the Euclidean distance between :math:`P(r, \varphi, \lambda)` and
+ :math:`Q(r_p, \varphi_p, \lambda_p)` can be calculated as follows [Grombein2013]_:
+
+ .. math::
+
+ l = \sqrt{ r^2 + r_p^2 - 2 r r_p \cos \Psi },
+
+ where
+
+ .. math::
+
+ \cos \Psi = \sin \varphi \sin \varphi_p +
+ \cos \varphi \cos \varphi_p \cos(\lambda - \lambda_p).
+
+ The radial component of the acceleration vector on a local North-oriented
+ system whose origin is located on the point :math:`P(r, \varphi, \lambda)`
+ is given by [Grombein2013]_:
+
+ .. math::
+
+ g_r(P) = \frac{G m}{l^3} (r_p \cos \Psi - r).
+
+ .. warning::
+
+ When working in Cartesian coordinates, the **z direction points downwards**,
+ i.e. positive and negative values represent points below and above the surface,
+ respectively.
+
 
  Parameters
  ----------
  coordinates : list or array
- List or array containing `longitude`, `latitude` and `radius` of computation
- points defined on a spherical geocentric coordinate system.
- Both `longitude` and `latitude` should be in degrees and `radius` in meters.
+ List or array containing the coordinates of computation points in the following
+ order: ``easting``, ``northing`` and ``down`` (if coordinates given in
+ Cartesian coordiantes), or ``longitude``, ``latitude`` and ``radius`` (if given
+ on a spherical geocentric coordinate system).
+ All ``easting``, ``northing`` and ``down`` should be in meters.
+ Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` in
+ meters.
  points : list or array
- List or array containing the coordinates of the point masses `longitude_p`,
- `latitude_p`, `radius_p` defined on a spherical geocentric coordinate system.
- Both `longitude_p` and `latitude_p` should be in degrees and `radius` in meters.
+ List or array containing the coordinates of the point masses in the following
+ order: ``easting``, ``northing`` and ``down`` (if coordinates given in
+ Cartesian coordiantes), or ``longitude``, ``latitude`` and ``radius`` (if given
+ on a spherical geocentric coordinate system).
+ All ``easting``, ``northing`` and ``down`` should be in meters.
+ Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` in
+ meters.
  masses : list or array
  List or array containing the mass of each point mass in kg.
  field : str
  Gravitational field that wants to be computed.
- The available fields are:
+ The available fields in Cartesian coordinates are:
+
+ - Gravitational potential: ``potential``
+ - Downward acceleration: ``g_z``
+
+ The available fields in spherical geocentric coordinates are:
 
  - Gravitational potential: ``potential``
  - Radial acceleration: ``g_r``
+
+ coordinate_system : str (optional)
+ Coordinate system of the coordinates of the computation points and the point
+ masses. Available coordinates systems: ``cartesian``, ``spherical``.
+ Default ``cartesian``.
  dtype : data-type (optional)
  Data type assigned to resulting gravitational field, and coordinates of point
  masses and computation points. Default to ``np.float64``.
@@ -37,50 +126,122 @@ def point_mass_gravity(coordinates, points, masses, field, dtype="float64"):
  Returns
  -------
  result : array
- Gravitational field generated by the `point_mass` on the computation points
- defined in `coordinates`.
+ Gravitational field generated by the ``point_mass`` on the computation points
+ defined in ``coordinates``.
  The potential is given in SI units, the accelerations in mGal and the Marussi
  tensor components in Eotvos.
  """
- kernels = {"potential": kernel_potential, "g_r": kernel_g_r}
- if field not in kernels:
+ # Organize dispatchers and kernel functions inside dictionaries
+ dispatchers = {
+ "cartesian": jit_point_mass_cartesian,
+ "spherical": jit_point_mass_spherical,
+ }
+ kernels = {
+ "cartesian": {"potential": kernel_potential_cartesian, "g_z": kernel_g_z},
+ "spherical": {"potential": kernel_potential_spherical, "g_r": kernel_g_r},
+ }
+ # Sanity checks for coordinate_system and field
+ if coordinate_system not in ("cartesian", "spherical"):
+ raise ValueError(
+ "Coordinate system {} not recognized".format(coordinate_system)
+ )
+ if field not in kernels[coordinate_system]:
  raise ValueError("Gravity field {} not recognized".format(field))
  # Figure out the shape and size of the output array
  cast = np.broadcast(*coordinates[:3])
  result = np.zeros(cast.size, dtype=dtype)
  # Prepare arrays to be passed to the jitted functions
- longitude, latitude, radius = (
- np.atleast_1d(i).ravel().astype(dtype) for i in coordinates[:3]
- )
- longitude_p, latitude_p, radius_p = (
- np.atleast_1d(i).ravel().astype(dtype) for i in points[:3]
- )
+ coordinates = (np.atleast_1d(i).ravel().astype(dtype) for i in coordinates[:3])
+ points = (np.atleast_1d(i).ravel().astype(dtype) for i in points[:3])
  masses = np.atleast_1d(masses).astype(dtype).ravel()
  # Compute gravitational field
- jit_point_mass_gravity(
- longitude,
- latitude,
- radius,
- longitude_p,
- latitude_p,
- radius_p,
- masses,
- result,
- kernels[field],
+ dispatchers[coordinate_system](
+ *coordinates, *points, masses, result, kernels[coordinate_system][field]
  )
  result *= GRAVITATIONAL_CONST
  # Convert to more convenient units
- if field == "g_r":
+ if field in ("g_r", "g_z"):
  result *= 1e5 # SI to mGal
  return result.reshape(cast.shape)
 
 
 @jit(nopython=True)
-def jit_point_mass_gravity(
+def jit_point_mass_cartesian(
+ easting, northing, down, easting_p, northing_p, down_p, masses, out, kernel
+): # pylint: disable=invalid-name
+ """
+ Compute gravity field of point masses on computation points in Cartesian coordinates
+
+ Parameters
+ ----------
+ easting, northing, down : 1d-arrays
+ Coordinates of computation points in Cartesian coordinate system.
+ easting_p, northing_p, down_p : 1d-arrays
+ Coordinates of point masses in Cartesian coordinate system.
+ masses : 1d-array
+ Mass of each point mass in SI units.
+ out : 1d-array
+ Array where the gravitational field on each computation point will be appended.
+ It must have the same size of ``easting``, ``northing`` and ``down``.
+ kernel : func
+ Kernel function that will be used to compute the gravity field on the
+ computation points.
+ """
+ for l in range(easting.size):
+ for m in range(easting_p.size):
+ out[l] += masses[m] * kernel(
+ easting[l], northing[l], down[l], easting_p[m], northing_p[m], down_p[m]
+ )
+
+
+@jit(nopython=True)
+def kernel_potential_cartesian(easting, northing, down, easting_p, northing_p, down_p):
+ """
+ Kernel function for potential gravity field in Cartesian coordinates
+ """
+ return 1 / _distance_cartesian(
+ [easting, northing, down], [easting_p, northing_p, down_p]
+ )
+
+
+@jit(nopython=True)
+def kernel_g_z(easting, northing, down, easting_p, northing_p, down_p):
+ """
+ Kernel function for downward component of gravity gradient in Cartesian coordinates
+ """
+ distance_sq = _distance_cartesian_sq(
+ [easting, northing, down], [easting_p, northing_p, down_p]
+ )
+ return (down_p - down) / distance_sq ** (3 / 2)
+
+
+@jit(nopython=True)
+def _distance_cartesian_sq(point_a, point_b):
+ """
+ Calculate the square distance between two points given in Cartesian coordinates
+ """
+ easting, northing, down = point_a[:]
+ easting_p, northing_p, down_p = point_b[:]
+ distance_sq = (
+ (easting - easting_p) ** 2 + (northing - northing_p) ** 2 + (down - down_p) ** 2
+ )
+ return distance_sq
+
+
+@jit(nopython=True)
+def _distance_cartesian(point_a, point_b):
+ """
+ Calculate the distance between two points given in Cartesian coordinates
+ """
+ return np.sqrt(_distance_cartesian_sq(point_a, point_b))
+
+
+@jit(nopython=True)
+def jit_point_mass_spherical(
  longitude, latitude, radius, longitude_p, latitude_p, radius_p, masses, out, kernel
 ): # pylint: disable=invalid-name
  """
- Compute gravity field of point masses on computation points.
+ Compute gravity field of point masses on computation points in spherical coordiantes
 
  Parameters
  ----------
@@ -123,11 +284,11 @@ def jit_point_mass_gravity(
 
 
 @jit(nopython=True)
-def kernel_potential(
+def kernel_potential_spherical(
  longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
  """
- Kernel function for potential gravity field
+ Kernel function for potential gravity field in spherical coordinates
  """
  coslambda = np.cos(longitude_p - longitude)
  cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda
@@ -140,7 +301,7 @@ def kernel_g_r(
  longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
 ):
  """
- Kernel function for radial component of gravity gradient
+ Kernel function for radial component of gravity gradient in spherical coordinates
  """
  coslambda = np.cos(longitude_p - longitude)
  cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda

harmonica/forward/tesseroid.py

@@ -6,7 +6,7 @@
 from numpy.polynomial.legendre import leggauss
 
 from ..constants import GRAVITATIONAL_CONST
-from .point_mass import jit_point_mass_gravity, kernel_potential, kernel_g_r
+from .point_mass import jit_point_mass_spherical, kernel_potential_spherical, kernel_g_r
 
 STACK_SIZE = 100
 MAX_DISCRETIZATIONS = 100000
@@ -110,7 +110,7 @@ def tesseroid_gravity(
  array(-112.54539933)
 
  """
- kernels = {"potential": kernel_potential, "g_r": kernel_g_r}
+ kernels = {"potential": kernel_potential_spherical, "g_r": kernel_g_r}
  if field not in kernels:
  raise ValueError("Gravity field {} not recognized".format(field))
  # Figure out the shape and size of the output array
@@ -258,7 +258,7 @@ def jit_tesseroid_gravity(
  weights,
  )
  # Compute gravity fields
- jit_point_mass_gravity(
+ jit_point_mass_spherical(
  computation_point[0:1],
  computation_point[1:2],
  computation_point[2:3],