Skip to content

Commit

Permalink
Move forward modeling examples to own category (#72)
Browse files Browse the repository at this point in the history
There will be many of these it might help to have them separated for
emphasis. Simplify the tesseroid example to only compute g_r to reduce
the run time of the example. Put on lower latitude to see the warping
effect of curvature.
  • Loading branch information
leouieda authored and santisoler committed Jul 18, 2019
1 parent 8222ed7 commit 301ca90
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 67 deletions.
2 changes: 2 additions & 0 deletions examples/forward/README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
Forward Modeling
----------------
43 changes: 43 additions & 0 deletions examples/forward/tesseroid.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
"""
Tesseroid Forward Modelling
===========================
Computing the gravitational fields generated by regional or global scale structures
require to take into account the curvature of the Earth. One common approach is to use
spherical prisms also known as tesseroids. We will compute the radial component of the
gravitational acceleration generated by a single tesseroid on a computation grid through
the :func:`harmonica.tesseroid_gravity` function.
"""
import harmonica as hm
import verde as vd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs


# Get default ellipsoid (WGS84) to obtain the mean Earth radius
ellipsoid = hm.get_ellipsoid()

# Define tesseroid with top surface at the mean Earth radius, thickness of 10km
# (bottom = top - thickness) and density of 2670kg/m^3
tesseroid = [-70, -50, -40, -20, ellipsoid.mean_radius - 10e3, ellipsoid.mean_radius]
density = 2670

# Define computation points on a regular grid at 100km above the mean Earth radius
coordinates = vd.grid_coordinates(
region=[-80, -40, -50, -10],
shape=(80, 80),
extra_coords=100e3 + ellipsoid.mean_radius,
)

# Compute the radial component of the acceleration
gravity = hm.tesseroid_gravity(coordinates, tesseroid, density, field="g_r")
print(gravity)

# Plot the gravitational field
fig = plt.figure(figsize=(8, 9))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree())
plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
ax.coastlines()
ax.set_title("Radial component of gravitational acceleration")
plt.show()
67 changes: 0 additions & 67 deletions examples/tesseroid.py

This file was deleted.

0 comments on commit 301ca90

Please sign in to comment.