diff --git a/doc/references.rst b/doc/references.rst index 15b82b5b9..ba43242bd 100644 --- a/doc/references.rst +++ b/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 `__ .. [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 `__ +.. [Blakely1995] Blakely, R. (1995). Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. doi:`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 `__ +.. [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 `__ .. [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 `__ .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press. diff --git a/examples/forward/point_mass.py b/examples/forward/point_mass.py new file mode 100644 index 000000000..4c8f91ab4 --- /dev/null +++ b/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() diff --git a/harmonica/forward/point_mass.py b/harmonica/forward/point_mass.py index f6a200e7c..5b12660cf 100644 --- a/harmonica/forward/point_mass.py +++ b/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 diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index 021d51f47..d3d679e53 100644 --- a/harmonica/forward/tesseroid.py +++ b/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], diff --git a/harmonica/tests/test_point_mass.py b/harmonica/tests/test_point_mass.py index 3aba53965..dfa1dd861 100644 --- a/harmonica/tests/test_point_mass.py +++ b/harmonica/tests/test_point_mass.py @@ -9,19 +9,106 @@ from ..forward.point_mass import point_mass_gravity -def test_invalid_field(): - "Check if an invalid gravitational field is passed as argument" +def test_invalid_coordinate_system(): + "Check if invalid coordinate system is passed" + coordinates = [0.0, 0.0, 0.0] point_mass = [0.0, 0.0, 0.0] mass = 1.0 - longitude = np.array(0.0) - latitude = np.array(0.0) - height = np.array(0.0) with pytest.raises(ValueError): point_mass_gravity( - [longitude, latitude, height], point_mass, mass, "this-field-does-not-exist" + coordinates, + point_mass, + mass, + "potential", + "this-is-not-a-valid-coordinate-system", ) +def test_invalid_field(): + "Check if an invalid gravitational field is passed as argument" + coordinates = [0.0, 0.0, 0.0] + point_mass = [0.0, 0.0, 0.0] + mass = 1.0 + for coordinate_system in ("spherical", "cartesian"): + with pytest.raises(ValueError): + point_mass_gravity( + coordinates, + point_mass, + mass, + "this-field-does-not-exist", + coordinate_system, + ) + + +def test_invalid_field_for_coordinate_system(): + "Check if an invalid field is passed for that coordinate system" + coordinates = [0.0, 0.0, 0.0] + point_mass = [0.0, 0.0, 0.0] + mass = 1.0 + cartesian_exclusive_fields = ["g_z"] + spherical_exclusive_fields = ["g_r"] + for field in spherical_exclusive_fields: + with pytest.raises(ValueError): + point_mass_gravity(coordinates, point_mass, mass, field, "cartesian") + for field in cartesian_exclusive_fields: + with pytest.raises(ValueError): + point_mass_gravity(coordinates, point_mass, mass, field, "spherical") + + +# --------------------------- +# Cartesian coordinates tests +# --------------------------- +@pytest.mark.use_numba +def test_potential_cartesian_symmetry(): + """ + Test if potential field of a point mass has symmetry in Cartesian coordinates + """ + # Define a single point mass + point_mass = [1.1, 1.2, 1.3] + masses = [2670] + # Define a set of computation points at a fixed distance from the point mass + distance = 3.3 + easting = point_mass[0] * np.ones(6) + northing = point_mass[1] * np.ones(6) + down = point_mass[2] * np.ones(6) + easting[0] += distance + easting[1] -= distance + northing[2] += distance + northing[3] -= distance + down[4] += distance + down[5] -= distance + coordinates = [easting, northing, down] + # Compute potential gravity field on each computation point + results = point_mass_gravity( + coordinates, point_mass, masses, "potential", "cartesian" + ) + npt.assert_allclose(*results) + + +@pytest.mark.use_numba +def test_g_z_symmetry(): + """ + Test if g_z field of a point mass has symmetry in Cartesian coordinates + """ + # Define a single point mass + point_mass = [1.1, 1.2, 1.3] + masses = [2670] + # Define a pair of computation points above and bellow the point mass + distance = 3.3 + easting = point_mass[0] * np.ones(2) + northing = point_mass[1] * np.ones(2) + down = point_mass[2] * np.ones(2) + down[0] += distance + down[1] -= distance + coordinates = [easting, northing, down] + # Compute g_z gravity field on each computation point + results = point_mass_gravity(coordinates, point_mass, masses, "g_z", "cartesian") + npt.assert_allclose(results[0], -results[1]) + + +# --------------------------- +# Spherical coordinates tests +# --------------------------- @pytest.mark.use_numba def test_point_mass_on_origin(): "Check potential and g_r of point mass on origin" @@ -39,7 +126,9 @@ def test_point_mass_on_origin(): # Compare results with analytical solutions for field in analytical: npt.assert_allclose( - point_mass_gravity([longitude, latitude, radius], point_mass, mass, field), + point_mass_gravity( + [longitude, latitude, radius], point_mass, mass, field, "spherical" + ), analytical[field], ) @@ -67,7 +156,9 @@ def test_point_mass_same_radial_direction(): # Compare results with analytical solutions for field in analytical: npt.assert_allclose( - point_mass_gravity(coordinates, point_mass, mass, field), + point_mass_gravity( + coordinates, point_mass, mass, field, "spherical" + ), analytical[field], ) @@ -95,7 +186,9 @@ def test_point_mass_potential_on_equator(): analytical = {"potential": GRAVITATIONAL_CONST * mass / distance} # Compare results with analytical solutions npt.assert_allclose( - point_mass_gravity(coordinates, point_mass, mass, "potential"), + point_mass_gravity( + coordinates, point_mass, mass, "potential", "spherical" + ), analytical["potential"], ) @@ -123,6 +216,8 @@ def test_point_mass_potential_on_same_meridian(): analytical = {"potential": GRAVITATIONAL_CONST * mass / distance} # Compare results with analytical solutions npt.assert_allclose( - point_mass_gravity(coordinates, point_mass, mass, "potential"), + point_mass_gravity( + coordinates, point_mass, mass, "potential", "spherical" + ), analytical["potential"], )