# Source code for geoana.em.fdem

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, epsilon_0
import numpy as np
import warnings
import properties

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

__all__ = [
'omega', 'wavenumber', 'skin_depth', 'sigma_hat', 'BaseFDEM',
'ElectricDipoleWholeSpace', 'MagneticDipoleWholeSpace'
]

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

[docs]def omega(frequency):
"""
Angular frequency

.. math::

\\omega = 2 \\pi f

**Required**
:param frequency float: frequency (Hz)

"""
return 2*np.pi*frequency

[docs]def wavenumber(
frequency, sigma, mu=mu_0, epsilon=epsilon_0, quasistatic=False
):
"""
Wavenumber of an electromagnetic wave in a medium with constant physical
properties

.. math::

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

**Required**

:param (float, numpy.ndarray) frequency: frequency (Hz)
: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
:param float epsilon: dielectric permittivity (F/m). Default: :math:\\epsilon_0 = 8.85 \\times 10^{-12} F/m
:param bool quasistatic: use the quasi-static assumption? Default: False

"""
w = omega(frequency)
if quasistatic is True:
return np.sqrt(-1j * w * mu * sigma)
return np.sqrt(w**2 * mu * epsilon - 1j * w * mu * sigma)

[docs]def skin_depth(frequency, sigma, mu=mu_0):
"""
Distance at which an em wave has decayed by a factor of :math:1/e in a
medium with constant physical properties

.. math::

\\sqrt{\\frac{2}{\\omega \\sigma \\mu}}

**Required**

:param float frequency: frequency (Hz)
: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

"""
w = omega(frequency)
return np.sqrt(2./(w*sigma*mu))

[docs]def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasistatic=False):
"""
conductivity with displacement current contribution

.. math::

\hat{\sigma} = \sigma + i \omega \\varepsilon

**Required**

:param (float, numpy.array) frequency: frequency (Hz)
:param float sigma: electrical conductivity (S/m)

**Optional**

:param float epsilon: dielectric permittivity. Default :math:\\varepsilon_0
:param bool quasistatic: use the quasi-static assumption? Default: False

"""
if quasistatic is True:
return sigma
return sigma + 1j*omega(frequency)*epsilon

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

[docs]class BaseFDEM(BaseEM):
"""
Base frequency domain electromagnetic class
"""
frequency = properties.Float(
"Source frequency (Hz)",
default=1.,
min=0.0
)

quasistatic = properties.Bool(
"Use the quasi-static approximation and ignore displacement current?",
default=False
)

@property
def omega(self):
"""
Angular frequency

.. math::

\\omega = 2\\pi f

"""
return omega(self.frequency)

@property
def sigma_hat(self):
"""
conductivity with displacement current contribution

.. math::

\\hat{\\sigma} = \\sigma + i \\omega \\varepsilon

"""
return sigma_hat(
self.frequency, self.sigma, epsilon=self.epsilon,
quasistatic=self.quasistatic
)

@property
def wavenumber(self):
"""
Wavenumber of an electromagnetic wave in a medium with constant
physical properties

.. math::

k = \\sqrt{\\omega**2 \\mu \\varepsilon - i \\omega \\mu \\sigma}

"""
return wavenumber(
self.frequency, self.sigma, mu=self.mu, epsilon=self.epsilon,
quasistatic=self.quasistatic
)

@property
def skin_depth(self):
"""
Distance at which an em wave has decayed by a factor of :math:1/e in
a medium with constant physical properties

.. math::

\\sqrt{\\frac{2}{\\omega \\sigma \\mu}}

"""
return skin_depth(self.frequency, self.sigma, mu=self.mu)

[docs]class ElectricDipoleWholeSpace(BaseElectricDipole, BaseFDEM):
"""
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 vector_potential(self, xyz):
"""
Vector potential for an electric dipole in a wholespace

.. math::

\\mathbf{A} = \\frac{I ds}{4 \\pi r} e^{-ikr}\\mathbf{\\hat{u}}

"""
r = self.distance(xyz)
a = (
(self.current * self.length) / (4*np.pi*r) *
np.exp(-i*self.wavenumber*r)
)
a = np.kron(np.ones(1, 3), np.atleast_2d(a).T)
return self.dot_orientation(a)

[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 = spatial.repeat_scalar(self.distance(xyz))
kr = self.wavenumber * r
ikr = 1j * kr

front_term = (
(self.current * self.length) / (4 * np.pi * self.sigma * r**3) *
np.exp(-ikr)
)
symmetric_term = (
spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz *
(-kr**2 + 3*ikr + 3) / r**2
)
oriented_term = (
(kr**2 - ikr - 1) *
np.kron(self.orientation, np.ones((dxyz.shape[0], 1)))
)
return front_term * (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

.. math::

\\mathbf{H} = \\nabla \\times \\mathbf{A}

"""
dxyz = self.vector_distance(xyz)
r = spatial.repeat_scalar(self.distance(xyz))
kr = self.wavenumber * r
ikr = 1j*kr

front_term = (
self.current * self.length / (4 * np.pi * r**2) * (ikr + 1) *
np.exp(-ikr)
)
return -front_term * self.cross_orientation(dxyz) / r

[docs]    def magnetic_flux_density(self, xyz):
"""
magnetic flux density from an electric dipole
"""
return self.mu * self.magnetic_field(xyz)

[docs]class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM):
"""
Harmonic magnetic dipole in a whole space.
"""

[docs]    def vector_potential(self, xyz):
"""
Vector potential for a magnetic dipole in a wholespace

.. math::

\\mathbf{F} = \\frac{i \\omega \\mu m}{4 \\pi r} e^{-ikr}
\\mathbf{\\hat{u}}

"""
r = self.distance(xyz)
f = (
(1j * self.omega * self.mu * self.moment) / (4 * np.pi * r) *
np.exp(-1j * self.wavenumber * r)
)
f = np.kron(np.ones(1, 3), np.atleast_2d(f).T)
return self.dot_orientation(f)

[docs]    def electric_field(self, xyz):
"""
Electric field from a magnetic dipole in a wholespace
"""
dxyz = self.vector_distance(xyz)
r = spatial.repeat_scalar(self.distance(xyz))
kr = self.wavenumber*r
ikr = 1j * kr

front_term = (
(1j * self.omega * self.mu * self.moment) / (4. * np.pi * r**2) *
(ikr + 1) * np.exp(-ikr)
)
return front_term * self.cross_orientation(dxyz) / r

[docs]    def current_density(self, xyz):
"""
Current density from a magnetic dipole in a wholespace
"""
return self.sigma * self.electric_field(xyz)

[docs]    def magnetic_field(self, xyz):
"""
Magnetic field due to a magnetic dipole in a wholespace
"""
dxyz = self.vector_distance(xyz)
r = spatial.repeat_scalar(self.distance(xyz))
kr = self.wavenumber*r
ikr = 1j*kr

front_term = self.moment / (4. * np.pi * r**3) * np.exp(-ikr)
symmetric_term = (
spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz *
(-kr**2 + 3*ikr + 3) / r**2
)
oriented_term = (
(kr**2 - ikr - 1) *
np.kron(self.orientation, np.ones((dxyz.shape[0], 1)))
)

return front_term * (symmetric_term + oriented_term)

[docs]    def magnetic_flux_density(self, xyz):
"""
Magnetic flux density due to a magnetic dipole in a wholespace
"""
return self.mu * self.magnetic_field(xyz)