geoana.em.fdem.ElectricDipoleWholeSpace.magnetic_field

ElectricDipoleWholeSpace.magnetic_field(xyz)

Magnetic field produced by the harmonic current dipole at a set of gridded locations.

For an electric current dipole oriented in the \(\hat{u}\) direction with dipole moment \(I ds\) and harmonic frequency \(f\), this method computes the magnetic field at the set of gridded xyz locations provided.

The analytic solution is adapted from Ward and Hohmann (1988). For a harmonic electric current dipole oriented in the \(\hat{x}\) direction, the solution at vector distance \(\mathbf{r}\) from the current dipole is:

\[\mathbf{H}(\mathbf{r}) = \frac{I ds}{4\pi r^2} (ikr + 1) e^{-ikr} \big ( - \frac{z}{r}\hat{y} + \frac{y}{r}\hat{z} \big )\]

where

\[k = \sqrt{\omega^2 \mu \varepsilon - i \omega \mu \sigma}\]

Parameters:

xyz(…, 3) numpy.ndarray

Gridded xyz locations

Returns:

(n_freq, …, 3) numpy.array of complex

Magnetic field at all frequencies for the gridded locations provided. Output array is squeezed when n_freq and/or n_loc = 1.

Examples

Here, we define an z-oriented electric dipole and plot the magnetic field on the xy-plane that intercepts z=0.

>>> from geoana.em.fdem import ElectricDipoleWholeSpace
>>> from geoana.utils import ndgrid
>>> from geoana.plotting_utils import plot2Ddata
>>> import numpy as np
>>> import matplotlib.pyplot as plt

Let us begin by defining the electric current dipole.

>>> frequency = np.logspace(1, 3, 3)
>>> location = np.r_[0., 0., 0.]
>>> orientation = np.r_[0., 0., 1.]
>>> current = 1.
>>> sigma = 1.0
>>> simulation = ElectricDipoleWholeSpace(
>>>     frequency, location=location, orientation=orientation,
>>>     current=current, sigma=sigma
>>> )

Now we create a set of gridded locations and compute the magnetic field.

>>> xyz = ndgrid(np.linspace(-1, 1, 20), np.linspace(-1, 1, 20), np.array([0]))
>>> H = simulation.magnetic_field(xyz)

Finally, we plot the real and imaginary components of the magnetic field.

>>> f_ind = 1
>>> fig = plt.figure(figsize=(6, 3))
>>> ax1 = fig.add_axes([0.15, 0.15, 0.40, 0.75])
>>> plot2Ddata(
>>>     xyz[:, 0:2], np.real(H[f_ind, :, 0:2]), vec=True, ax=ax1, scale='log', ncontour=25
>>> )
>>> ax1.set_xlabel('X')
>>> ax1.set_ylabel('Y')
>>> ax1.autoscale(tight=True)
>>> ax1.set_title('Real component {} Hz'.format(frequency[f_ind]))
>>> ax2 = fig.add_axes([0.6, 0.15, 0.40, 0.75])
>>> plot2Ddata(
>>>     xyz[:, 0:2], np.imag(H[f_ind, :, 0:2]), vec=True, ax=ax2, scale='log', ncontour=25
>>> )
>>> ax2.set_xlabel('X')
>>> ax2.set_yticks([])
>>> ax2.autoscale(tight=True)
>>> ax2.set_title('Imag component {} Hz'.format(frequency[f_ind]))

(Source code, png, pdf)