Magnetic dipoles#
We can compute the magnetic field generated by dipole sources through
the harmonica.dipole_magnetic
function.
Each dipole can be defined as a tuple containing its coordinates in the following order: easting, northing, upward (in Cartesian coordinates and in meters).
Let’s define a single dipole and compute the magnetic field it generates on a computation point above it.
Define the dipole and its magnetic moment vector:
dipole = (20, 40, -50)
magnetic_moment = (100, 100, 100)
Define the observation point above it:
coordinates = (20, 40, 10)
And compute the three components of the magnetic field the dipole generates by
choosing field="b"
:
b_e, b_n, b_u = hm.dipole_magnetic(coordinates, dipole, magnetic_moment, field="b")
print(b_e, b_n, b_u)
-0.04629629632149888 -0.04629629632149888 0.09259259264299775
We can compute just a single component with field
set to "b_e"
,
"b_n"
or "b_u"
:
b_u = hm.dipole_magnetic(
coordinates, dipole, magnetic_moment, field="b_u"
)
print(b_u)
0.09259259264299775
Multiple dipoles#
We can also use this function to compute the magnetic fields of multiple dipoles on multiple observation points.
Let’s define a regular grid of observation points:
import verde as vd
region = (-100, 100, -100, 100)
spacing = 1
height = 0
coordinates = vd.grid_coordinates(
region=region, spacing=spacing, extra_coords=height
)
And a set of dipoles with their magnetic moments:
import numpy as np
easting = [25, 35, -30, -50]
northing = [3, -38, 22, -30]
upward = [-200, -100, -300, -150]
dipoles = (easting, northing, upward)
mag_e = [1e3, 2e3, 500, 2e3]
mag_n = [1e3, 2e3, 500, 2e3]
mag_u = [1e3, 2e3, 500, 2e3]
magnetic_moments = (mag_e, mag_n, mag_u)
Now, let’s compute the magnetic field components that the dipoles generate on every observation point:
b_e, b_n, b_u = hm.dipole_magnetic(coordinates, dipoles, magnetic_moments, field="b")
import matplotlib.pyplot as plt
fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
fields = {"b_e": b_e, "b_n": b_n, "b_u": b_u}
for field, ax in zip(fields, axes):
tmp = ax.pcolormesh(coordinates[0], coordinates[1], fields[field])
ax.set_aspect("equal")
ax.set_title(field)
ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="both")
plt.colorbar(tmp, ax=ax, orientation="horizontal", label="nT", pad=0.008)
plt.show()