# geoana.em.fdem.ElectricDipoleWholeSpace.current_density#

ElectricDipoleWholeSpace.current_density(xyz)#

Current density for 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 current density 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:

$\begin{split}\;\;\;\; \mathbf{J}(\mathbf{r}) = \frac{\sigma I ds}{4\pi (\sigma + i \omega \varepsilon) r^3} & e^{-ikr} ... \\ & \Bigg [ \Bigg ( \frac{x^2}{r^2}\hat{x} + \frac{xy}{r^2}\hat{y} + \frac{xz}{r^2}\hat{z} \Bigg ) \big ( -k^2r^2 + 3ikr + 3 \big ) \big ( k^2 r^2 - ikr - 1 \big ) \hat{x} \Bigg ]\end{split}$

where

$k = \sqrt{\omega^2 \mu \varepsilon - i \omega \mu \sigma}$
Parameters:
xyz(n, 3) numpy.ndarray

Gridded xyz locations

Returns:
(n_freq, n_loc, 3) numpy.array of complex

Current density 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 x-oriented electric dipole and plot the current density on the xz-plane that intercepts y=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_[1., 0., 0.]
>>> 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 current density.

>>> xyz = ndgrid(np.linspace(-1, 1, 20), np.array(), np.linspace(-1, 1, 20))
>>> J = simulation.current_density(xyz)


Finally, we plot the real and imaginary components of the current density.

>>> 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(J[f_ind, :, 0::2]), vec=True, ax=ax1, scale='log', ncontour=25
>>> )
>>> ax1.set_xlabel('X')
>>> ax1.set_ylabel('Z')
>>> 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(J[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]))