|
| 1 | +""" |
| 2 | +Forward modelling for point masses |
| 3 | +""" |
| 4 | +import numpy as np |
| 5 | +from numba import jit |
| 6 | + |
| 7 | +from ..constants import GRAVITATIONAL_CONST |
| 8 | + |
| 9 | + |
| 10 | +def point_mass_gravity(coordinates, points, masses, field, dtype="float64"): |
| 11 | + """ |
| 12 | + Compute gravitational fields of a point mass on spherical coordinates. |
| 13 | +
|
| 14 | + Parameters |
| 15 | + ---------- |
| 16 | + coordinates : list or array |
| 17 | + List or array containing `longitude`, `latitude` and `radius` of computation |
| 18 | + points defined on a spherical geocentric coordinate system. |
| 19 | + Both `longitude` and `latitude` should be in degrees and `radius` in meters. |
| 20 | + points : list or array |
| 21 | + List or array containing the coordinates of the point masses `longitude_p`, |
| 22 | + `latitude_p`, `radius_p` defined on a spherical geocentric coordinate system. |
| 23 | + Both `longitude_p` and `latitude_p` should be in degrees and `radius` in meters. |
| 24 | + masses : list or array |
| 25 | + List or array containing the mass of each point mass in kg. |
| 26 | + field : str |
| 27 | + Gravitational field that wants to be computed. |
| 28 | + The available fields are: |
| 29 | +
|
| 30 | + - Gravitational potential: ``potential`` |
| 31 | + - Radial acceleration: ``g_r`` |
| 32 | + dtype : data-type (optional) |
| 33 | + Data type assigned to resulting gravitational field, and coordinates of point |
| 34 | + masses and computation points. Default to ``np.float64``. |
| 35 | +
|
| 36 | +
|
| 37 | + Returns |
| 38 | + ------- |
| 39 | + result : array |
| 40 | + Gravitational field generated by the `point_mass` on the computation points |
| 41 | + defined in `coordinates`. |
| 42 | + The potential is given in SI units, the accelerations in mGal and the Marussi |
| 43 | + tensor components in Eotvos. |
| 44 | + """ |
| 45 | + kernels = {"potential": kernel_potential, "g_r": kernel_g_r} |
| 46 | + if field not in kernels: |
| 47 | + raise ValueError("Gravity field {} not recognized".format(field)) |
| 48 | + # Figure out the shape and size of the output array |
| 49 | + cast = np.broadcast(*coordinates[:3]) |
| 50 | + result = np.zeros(cast.size, dtype=dtype) |
| 51 | + # Prepare arrays to be passed to the jitted functions |
| 52 | + longitude, latitude, radius = ( |
| 53 | + np.atleast_1d(i).ravel().astype(dtype) for i in coordinates[:3] |
| 54 | + ) |
| 55 | + longitude_p, latitude_p, radius_p = ( |
| 56 | + np.atleast_1d(i).ravel().astype(dtype) for i in points[:3] |
| 57 | + ) |
| 58 | + masses = np.atleast_1d(masses).astype(dtype).ravel() |
| 59 | + # Compute gravitational field |
| 60 | + jit_point_mass_gravity( |
| 61 | + longitude, |
| 62 | + latitude, |
| 63 | + radius, |
| 64 | + longitude_p, |
| 65 | + latitude_p, |
| 66 | + radius_p, |
| 67 | + masses, |
| 68 | + result, |
| 69 | + kernels[field], |
| 70 | + ) |
| 71 | + result *= GRAVITATIONAL_CONST |
| 72 | + # Convert to more convenient units |
| 73 | + if field == "g_r": |
| 74 | + result *= 1e5 # SI to mGal |
| 75 | + return result.reshape(cast.shape) |
| 76 | + |
| 77 | + |
| 78 | +@jit(nopython=True) |
| 79 | +def jit_point_mass_gravity( |
| 80 | + longitude, latitude, radius, longitude_p, latitude_p, radius_p, masses, out, kernel |
| 81 | +): # pylint: disable=invalid-name |
| 82 | + """ |
| 83 | + Compute gravity field of point masses on computation points. |
| 84 | +
|
| 85 | + Parameters |
| 86 | + ---------- |
| 87 | + longitude, latitude, radius : 1d-arrays |
| 88 | + Coordinates of computation points in spherical geocentric coordinate system. |
| 89 | + longitude_p, latitude_p, radius_p : 1d-arrays |
| 90 | + Coordinates of point masses in spherical geocentric coordinate system. |
| 91 | + masses : 1d-array |
| 92 | + Mass of each point mass in SI units. |
| 93 | + out : 1d-array |
| 94 | + Array where the gravitational field on each computation point will be appended. |
| 95 | + It must have the same size of ``longitude``, ``latitude`` and ``radius``. |
| 96 | + kernel : func |
| 97 | + Kernel function that will be used to compute the gravity field on the |
| 98 | + computation points. |
| 99 | + """ |
| 100 | + # Compute quantities related to computation point |
| 101 | + longitude = np.radians(longitude) |
| 102 | + latitude = np.radians(latitude) |
| 103 | + cosphi = np.cos(latitude) |
| 104 | + sinphi = np.sin(latitude) |
| 105 | + # Compute quantities related to point masses |
| 106 | + longitude_p = np.radians(longitude_p) |
| 107 | + latitude_p = np.radians(latitude_p) |
| 108 | + cosphi_p = np.cos(latitude_p) |
| 109 | + sinphi_p = np.sin(latitude_p) |
| 110 | + # Compute gravity field |
| 111 | + for l in range(longitude.size): |
| 112 | + for m in range(longitude_p.size): |
| 113 | + out[l] += masses[m] * kernel( |
| 114 | + longitude[l], |
| 115 | + cosphi[l], |
| 116 | + sinphi[l], |
| 117 | + radius[l], |
| 118 | + longitude_p[m], |
| 119 | + cosphi_p[m], |
| 120 | + sinphi_p[m], |
| 121 | + radius_p[m], |
| 122 | + ) |
| 123 | + |
| 124 | + |
| 125 | +@jit(nopython=True) |
| 126 | +def kernel_potential( |
| 127 | + longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p |
| 128 | +): |
| 129 | + """ |
| 130 | + Kernel function for potential gravity field |
| 131 | + """ |
| 132 | + coslambda = np.cos(longitude_p - longitude) |
| 133 | + cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda |
| 134 | + distance_sq = (radius - radius_p) ** 2 + 2 * radius * radius_p * (1 - cospsi) |
| 135 | + return 1 / np.sqrt(distance_sq) |
| 136 | + |
| 137 | + |
| 138 | +@jit(nopython=True) |
| 139 | +def kernel_g_r( |
| 140 | + longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p |
| 141 | +): |
| 142 | + """ |
| 143 | + Kernel function for radial component of gravity gradient |
| 144 | + """ |
| 145 | + coslambda = np.cos(longitude_p - longitude) |
| 146 | + cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda |
| 147 | + distance_sq = (radius - radius_p) ** 2 + 2 * radius * radius_p * (1 - cospsi) |
| 148 | + delta_z = radius_p * cospsi - radius |
| 149 | + return delta_z / distance_sq ** (3 / 2) |
0 commit comments