# Rectangular prisms

## Contents

# Rectangular prisms¶

One of the geometric bodies that Harmonica is able to forward model are
*rectangular prisms*.
We can compute their gravitational fields through the
`harmonica.prism_gravity`

function.

Rectangular prisms can only be defined in Cartesian coordinates, therefore they are intended to be used to forward model bodies in a small geographic region, in which we can neglect the curvature of the Earth.

Each prism can be defined as a tuple containing its boundaries in the following
order: *west*, *east*, *south*, *north*, *bottom*, *top* in Cartesian
coordinates and in meters.

Lets define a single prism and compute the gravitational potential it generates on a computation point located at 500 meters above its uppermost surface.

Define the prism and its density (in kg per cubic meters):

```
prism = (-100, 100, -100, 100, -1000, 500)
density = 3000
```

Define a computation point located 500 meters above the *top* surface of the
prism:

```
coordinates = (0, 0, 1000)
```

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

```
potential = hm.prism_gravity(coordinates, prism, density, field="potential")
print(potential, "J/kg")
```

```
0.011051569953086739 J/kg
```

## Passing multiple prisms¶

We can compute the gravitational field of a set of prisms by passing a list of
them, where each prism is defined as mentioned before, and then making
a single call of the `harmonica.prism_gravity`

function.

Lets define a set of four prisms, along with their densities:

```
prisms = [
[2e3, 3e3, 2e3, 3e3, -10e3, -1e3],
[3e3, 4e3, 7e3, 8e3, -9e3, -1e3],
[7e3, 8e3, 1e3, 2e3, -7e3, -1e3],
[8e3, 9e3, 6e3, 7e3, -8e3, -1e3],
]
densities = [2670, 3300, 2900, 2980]
```

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=(0, 10e3, 0, 10e3), shape=(40, 40), extra_coords=0
)
```

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

```
g_z = hm.prism_gravity(coordinates, prisms, densities, field="g_z")
```

Note

When passing multiple prisms and coordinates to
`harmonica.prism_gravity`

we calculate the field in parallel using
multiple CPUs, speeding up the computation.

Lets plot this gravitational field:

```
import pygmt
grid = vd.make_xarray_grid(
coordinates, g_z, data_names="g_z", extra_coords_names="extra")
fig = pygmt.Figure()
fig.grdimage(
region=(0, 10e3, 0, 10e3),
projection="X10c",
grid=grid.g_z,
frame=["WSne", "x+leasting (m)", "y+lnorthing (m)"],
cmap='viridis',)
fig.colorbar(cmap=True, position="JMR", frame=["a2", "x+lmGal"])
fig.show()
```

## Prism layer¶

One of the most common usage of prisms is to model geologic structures.
Harmonica offers the possibility to define a layer of prisms through the
`harmonica.prism_layer`

function: a regular grid of
prisms of equal size along the horizontal dimensions and with variable top and
bottom boundaries.
It returns a `xarray.Dataset`

with the coordinates of the centers of the
prisms and their corresponding physical properties.

The `harmonica.DatasetAccessorPrismsLayer`

Dataset accessor can be used
to obtain some properties of the layer like its shape and size or the
boundaries of any prism in the layer.
Moreover, we can use the `harmonica.DatasetAcessorPrismsLayer.gravity`

method to compute the gravitational field of the prism layer on any set of
computation points.

Lets create a simple prism layer, whose lowermost boundaries will be set on zero and their uppermost boundary will approximate a sinusoidal function. We can start by setting the region of the layer and the horizontal dimensions of the prisms:

```
region = (0, 100e3, -40e3, 40e3)
spacing = 2000
```

Then we can define a regular grid where the centers of the prisms will fall:

```
easting, northing = vd.grid_coordinates(region=region, spacing=spacing)
```

We need to define a 2D array for the uppermost *surface* of the layer. We will
sample a trigonometric function for this simple example:

```
import numpy as np
wavelength = 24 * spacing
surface = np.abs(np.sin(easting * 2 * np.pi / wavelength))
```

Lets assign the same density to each prism through a 2d array with the same value: 2700 kg per cubic meter.

```
density = np.full_like(surface, 2700)
```

Now we can define the prism layer specifying the reference level to zero:

```
prisms = hm.prism_layer(
coordinates=(easting, northing),
surface=surface,
reference=0,
properties={"density": density},
)
```

Lets define a grid of observation points at 1 km above the zeroth height:

```
region_pad = vd.pad_region(region, 10e3)
coordinates = vd.grid_coordinates(
region_pad, spacing=spacing, extra_coords=1e3
)
```

And compute the gravitational field generated by the prism layer on them:

```
gravity = prisms.prism_layer.gravity(coordinates, field="g_z")
```

Finally, lets plot the gravitational field:

```
grid = vd.make_xarray_grid(
coordinates, gravity, data_names="gravity", extra_coords_names="extra")
fig = pygmt.Figure()
title = "Gravitational acceleration of a layer of prisms"
fig.grdimage(
region=region_pad,
projection="X10c",
grid=grid.gravity,
frame=[f"WSne+t{title}", "x+leasting (m)", "y+lnorthing (m)"],
cmap='viridis',)
fig.colorbar(cmap=True, position="JMR", frame=["a.02", "x+lmGal"])
fig.show()
```