# Point sources¶

We can compute the gravitational field of point masses through the `harmonica.point_gravity` function. It offers the possibility to define point sources and computation points either in Cartesian or in spherical coordinates.

Each point mass can be defined as a tuple containing its coordinates in the following order: easting, northing and upward (in Cartesian coordinates) or longitude, spherical latitude and radius (in spherical coordinates).

## Cartesian coordinates¶

For example, lets define a single point mass and compute the gravitational potential it generates on a computation point located 100 meters above it.

Define the point source and its mass (in kg):

```point = (400, 300, 200)
mass = 1e6
```

Define a computation point located 100m above the point source:

```coordinates = (point, point, point + 100)
```

And finally, compute the gravitational potential field it generates on the computation point:

```potential = hm.point_gravity(coordinates, point, mass, field="potential")
print(potential, "J/kg")
```
```6.672999999999999e-07 J/kg
```

### Working with multiple sources¶

We can easily compute the gravitational field of a set of point sources by grouping its coordinates in arrays and with a single call of the `harmonica.point_gravity`.

Lets define a set of four point sources, along with their masses:

```import numpy as np

easting = np.array([250, 750, 250, 750])
northing = np.array([250, 250, 750, 750])
upward = np.array([-100, -100, -100, -100])
points = (easting, northing, upward)

masses = np.array([1e6, -1e6, 2e6, -3e6])
```

We can define a set of computation points located on a regular grid at zero height:

```import verde as vd

coordinates = vd.grid_coordinates(
region=(-250, 1250, -250, 1250), shape=(40, 40), extra_coords=0
)
```

And finally calculate the vertical component of the gravitational acceleration generated by the whole set of point sources on every computation point:

```g_z = hm.point_gravity(coordinates, points, masses, field="g_z")
```

Note

When passing multiple sources and coordinates to `harmonica.point_gravity` we calculate the field in parallel using multiple CPUs, speeding up the computation.

Lets plot this gravitational field:

```import matplotlib.pyplot as plt

maxabs = vd.maxabs(g_z)
plt.pcolormesh(
*coordinates[:2], g_z, vmin=-maxabs, vmax=maxabs, cmap="seismic"
)
plt.colorbar(label="mGal")
plt.gca().set_aspect("equal")
plt.xlabel("easting [m]")
plt.ylabel("northing [m]")
plt.show()
``` ## Spherical coordinates¶

Alternatively, we can compute the gravitational fields of point sources defined in a spherical coordinate system. To do so, we need to pass the `coordinate_system` argument as `"spherical"``. The coordinates of the point source must now be passed as longitude, spherical latitude and radius, where the two former ones must be in decimal degrees and the latter in meters.

Lets define a single point source in the equator, at a longitude of 45 degrees and located on the surface of the WGS84 ellipsoid.

We can start by obtaining the WGS84 ellipsoid from `boule`:

```import boule as bl

ellipsoid = bl.WGS84
```

Then we can define the point source in the equator along with its mass:

```longitude, latitude = 45, 0
point = (longitude, latitude, radius)

mass = 1e6
```

We can also define a computation point located 1km above the point source:

```coordinates = (point, point, point + 1000)
```

And finally compute the gravitational field generated by the source on the computation point:

```g_z = hm.point_gravity(
coordinates, point, mass, field="g_z", coordinate_system="spherical"
)
print(g_z, "mGal")
```
```6.6729999999999995e-06 mGal
```

### Working with sources defined in geodetic coordinates¶

If our point sources and computation points are defined in geodetic coordinates, we can use the `boule.Ellipsoid.geodetic_to_spherical` method to convert them to spherical coordinates and use them to compute their gravitational field.

Lets define a set point sources in geodetic coordinates and their masses:

```longitude = np.array([-71, -71, -69, -69])
latitude = np.array([-45, -43, -45, -43])
height = np.array([-10e3, -20e3, -30e3, -20e3])
points = (longitude, latitude, height)

masses = np.array([1e6, 2e6, -3e6, 5e6])
```

Then, build a regular grid of computation points defined in geodetic coordinates and located at 1km above the ellipsoid:

```coordinates = vd.grid_coordinates(
region=(-72, -68, -46, -42),
shape=(101, 101),
extra_coords=20e3,
)
```

Before we can start forward modelling these sources, we need to convert them to spherical coordinates. To do so, we can use the `boule.Ellipsoid.geodetic_to_spherical` method:

```points_spherical = ellipsoid.geodetic_to_spherical(*points)
coordinates_spherical = ellipsoid.geodetic_to_spherical(*coordinates)
```

We can finally use these converted coordinates to compute the gravitational field the source generate on every computation point:

```g_z = hm.point_gravity(
coordinates_spherical,
points_spherical,
masses,
field="g_z",
coordinate_system="spherical",
)
```

Lets plot these results using `cartopy`:

```import cartopy.crs as ccrs

plt.figure(figsize=(8, 6))
ax = plt.axes(projection=ccrs.Mercator())
maxabs = vd.maxabs(g_z)
tmp = ax.pcolormesh(
*coordinates[:2],
g_z,
vmin=-maxabs,
vmax=maxabs,
cmap="seismic",
transform=ccrs.PlateCarree(),
)
ax.set_title("Gravitational acceleration (downward)")
plt.colorbar(tmp, ax=ax, pad=0.04, shrink=0.73, label="mGal")
plt.show()
``` 