# geoana.em.tdem.ElectricDipoleWholeSpace.current_density#

ElectricDipoleWholeSpace.current_density(xyz)#

Current density for the transient current dipole at a set of gridded locations.

For an electric current dipole oriented in the $$\hat{u}$$ direction with dipole moment $$I ds$$, this method computes the transient current density at the set of observation times for the gridded xyz locations provided.

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

$\begin{split}{\bf e_e}(t) = \frac{Ids}{4\pi \sigma r^2} & \Bigg [ \Bigg ( \frac{x^2}{r^2}\mathbf{\hat x} + \frac{xy}{r^2}\mathbf{\hat y} + \frac{xz}{r^2}\mathbf{\hat z} \Bigg ) \Bigg ( 3 \, \textrm{erf}(\theta r) - \bigg ( \frac{4}{\sqrt{\pi}}\theta^3 r^3 \; ... \\ & + \frac{6}{\sqrt{\pi}} \theta r \bigg ) e^{-\theta^2 r^2} \Bigg ) - \Bigg ( \textrm{erf}(\theta r) - \bigg ( \frac{4}{\sqrt{\pi}} \theta^3 r^3 + \frac{2}{\sqrt{\pi}} \theta r \bigg ) e^{-\theta^2 r^2} \Bigg ) \mathbf{\hat x} \Bigg ]\end{split}$

where

$\theta = \Bigg ( \frac{\mu\sigma}{4t} \Bigg )^{1/2}$
Parameters:
xyz(n, 3) numpy.ndarray

Gridded xyz locations

Returns:
(n_time, n_loc, 3) numpy.array of float

Current density at all times for the gridded locations provided. Output array is squeezed when n_time 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.tdem 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.

>>> time = np.logspace(-6, -2, 3)
>>> location = np.r_[0., 0., 0.]
>>> orientation = np.r_[1., 0., 0.]
>>> current = 1.
>>> sigma = 1.0
>>> simulation = ElectricDipoleWholeSpace(
>>>     time, 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(-10, 10, 20), np.array(), np.linspace(-10, 10, 20))
>>> J = simulation.current_density(xyz)


Finally, we plot the current density at the desired locations/times.

>>> t_ind = 0
>>> fig = plt.figure(figsize=(4, 4))
>>> ax = fig.add_axes([0.15, 0.15, 0.8, 0.8])
>>> plot2Ddata(xyz[:, 0::2], J[t_ind, :, 0::2], ax=ax, vec=True, scale='log')
>>> ax.set_xlabel('X')
>>> ax.set_ylabel('Z')
>>> ax.set_title('Current density at {} s'.format(time[t_ind]))