# Source code for geoana.em.tdem

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals

from scipy.constants import mu_0
from scipy.special import erf
import numpy as np
import properties

from .base import BaseElectricDipole, BaseEM
from .. import spatial

###############################################################################
#                                                                             #
#                           Utility Functions                                 #
#                                                                             #
###############################################################################

[docs]def peak_time(z, sigma, mu=mu_0):
"""
Peak time <https://em.geosci.xyz/content/maxwell1_fundamentals/transient_planewaves_homogeneous/peaktime.html>_:
Time at which the maximum signal amplitude is observed at a particular
location for a transient plane wave through a homogeneous medium.

**Required**

:param float z: distance from source (m)
:param float sigma: electrical conductivity (S/m)

**Optional**

:param float mu: magnetic permeability (H/m). Default: :math:\mu_0 = 4\pi \\times 10^{-7} H/m

"""
return (mu * sigma * z**2)/6.

[docs]def diffusion_distance(time, sigma, mu=mu_0):
"""
Diffusion distance <https://em.geosci.xyz/content/maxwell1_fundamentals/transient_planewaves_homogeneous/peakdistance.html>_:
Distance at which the signal amplitude is largest for a given time after
shut off. Also referred to as the peak distance
"""
return np.sqrt(2*time/(mu*sigma))

[docs]def theta(time, sigma, mu=mu_0):
"""
Analog to wavenumber in the frequency domain. See Ward and Hohmann, 1988
pages 174-175
"""
return np.sqrt(mu*sigma/(4.*time))

###############################################################################
#                                                                             #
#                                  Classes                                    #
#                                                                             #
###############################################################################

[docs]class BaseTDEM(BaseEM):

time = properties.Float(
"time after shut-off at which we are evaluating the fields (s)",
required=True,
default=1e-4
)

[docs]    def peak_time(self, z):
return peak_time(z, self.sigma, self.mu)

@property
def diffusion_distance(self):
return diffusion_distance(self.time, self.sigma, self.mu)

@property
def theta(self):
return theta(self.time, self.sigma, self.mu)

[docs]class ElectricDipoleWholeSpace(BaseElectricDipole, BaseTDEM):
"""
Harmonic electric dipole in a whole space. The source is
(c.f. Ward and Hohmann, 1988 page 173). The source current
density for a dipole located at :math:\mathbf{r}_s with orientation
:math:\mathbf{\hat{u}}

.. math::

\mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r}
- \mathbf{r}_s)\mathbf{\hat{u}}

"""

[docs]    def electric_field(self, xyz):
"""
Electric field from an electric dipole

.. math::

\mathbf{E} = \\frac{1}{\hat{\sigma}} \\nabla \\nabla \cdot \mathbf{A}
- i \omega \mu \mathbf{A}

"""
dxyz = self.vector_distance(xyz)
r = self.distance(xyz)
r = spatial.repeat_scalar(r)
thetar = self.theta * r
root_pi = np.sqrt(np.pi)

front = (
(self.current * self.length) / (4 * np.pi * self.sigma * r**3)
)

symmetric_term = (
(
- (
4/root_pi * thetar ** 3 + 6/root_pi * thetar
) * np.exp(-thetar**2) +
3 * erf(thetar)
) * (
spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz / r**2
)
)

oriented_term = (
(
4./root_pi * thetar**3 + 2./root_pi * thetar
) * np.exp(-thetar**2) -
erf(thetar)
) * np.kron(self.orientation, np.ones((dxyz.shape[0], 1)))

return front * (symmetric_term + oriented_term)

[docs]    def current_density(self, xyz):
"""
Current density due to a harmonic electric dipole
"""
return self.sigma * self.electric_field(xyz)

[docs]    def magnetic_field(self, xyz):
"""
Magnetic field from an electric dipole
"""
dxyz = self.vector_distance(xyz)
r = self.distance(dxyz)
r = spatial.repeat_scalar(r)
thetar = self.theta * r

front = (
self.current * self.length / (4 * np.pi * r**2) * (
2/root_pi * thetar * np.exp(-thetar**2) + erf(thetar)
)
)

return - front * self.cross_orientation(xyz) / r

[docs]    def magnetic_field_time_deriv(self, xyz):
"""
Time derivative of the magnetic field,
:math:\\frac{\partial \mathbf{h}}{\partial t}
"""

dxyz = self.vector_distance(xyz)
r = self.distance(xyz)
r = spatial.repeat_scalar(r)

front = (
self.current * self.length * self.theta**3 * r /
(2 * np.sqrt(np.pi)**3 * self.time)
)

return - front * self.cross_orientation(xyz) / r

[docs]    def magnetic_flux_density(self, xyz):
"""
Magnetic flux density from an electric dipole
"""

return self.mu * self.magnetic_field(xyz)

[docs]    def magnetic_flux_density_time_deriv(self, xyz):
"""
Time derivative of the magnetic flux density from an electric dipole
"""

return self.mu * self.magnetic_field_time_deriv(xyz)