from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
import numpy as np
import warnings
import properties
from scipy.special import ellipk, ellipe
from scipy.constants import mu_0
from .base import BaseEM, BaseDipole, BaseMagneticDipole, BaseElectricDipole
from .. import spatial
__all__ = ["MagneticDipoleWholeSpace", "CircularLoopWholeSpace"]
[docs]class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseEM):
"""
Static magnetic dipole in a wholespace.
"""
[docs] def vector_potential(self, xyz, coordinates="cartesian"):
"""Vector potential of a static magnetic dipole. See Griffiths, 1999
equation 5.83
.. math::
\\vec{A}(\\vec{r}) = \\frac{\mu_0}{4\pi}
\\frac{\\vec{m}\\times\\vec{r}}{r^3}
**Required**
:param numpy.ndarray xyz: Location at which we calculate the vector
potential
**Optional**
:param str coordinates: coordinate system that the xyz is provided
in and that the solution will be returned
in (cartesian or cylindrical).
Default: `"cartesian"`
**Returns**
:rtype: numpy.ndarray
:return: The magnetic vector potential at each observation location
"""
supported_coordinates = ["cartesian", "cylindrical"]
assert coordinates.lower() in supported_coordinates, (
"coordinates must be in {}, the coordinate system "
"you provided, {}, is not yet supported".format(
supported_coordinates, coordinates
)
)
n_obs = xyz.shape[0]
# orientation of the dipole
if coordinates.lower() == "cylindrical":
xyz = spatial.cylindrical_2_cartesian(xyz)
dxyz = self.vector_distance(xyz)
r = spatial.repeat_scalar(self.distance(xyz))
m = self.moment * np.atleast_2d(self.orientation).repeat(n_obs, axis=0)
m_cross_r = np.cross(m, dxyz)
a = (self.mu / (4 * np.pi)) * m_cross_r / (r**3)
if coordinates.lower() == "cylindrical":
a = spatial.cartesian_2_cylindrical(xyz, a)
return a
[docs] def magnetic_flux_density(self, xyz, coordinates="cartesian"):
"""Magnetic flux (:math:`\\vec{b}`) of a static magnetic dipole
**Required**
:param numpy.ndarray xyz: Location of the receivers(s)
**Optional**
:param str coordinates: coordinate system that the xyz is provided
in and that the solution will be returned
in (cartesian or cylindrical).
Default: `"cartesian"`
**Returns**
:rtype: numpy.ndarray
:return: The magnetic flux at each observation location
"""
supported_coordinates = ["cartesian", "cylindrical"]
assert coordinates.lower() in supported_coordinates, (
"coordinates must be in {}, the coordinate system "
"you provided, {}, is not yet supported".format(
supported_coordinates, coordinates
)
)
n_obs = xyz.shape[0]
if coordinates.lower() == "cylindrical":
xyz = spatial.cylindrical_2_cartesian(xyz)
r = self.vector_distance(xyz)
dxyz = spatial.repeat_scalar(self.distance(xyz))
m_vec = (
self.moment * np.atleast_2d(self.orientation).repeat(n_obs, axis=0)
)
m_dot_r = (m_vec * r).sum(axis=1)
# Repeat the scalars
m_dot_r = np.atleast_2d(m_dot_r).T.repeat(3, axis=1)
# dxyz = np.atleast_2d(dxyz).T.repeat(3, axis=1)
b = (self.mu / (4 * np.pi)) * (
(3.0 * r * m_dot_r / (dxyz ** 5)) -
m_vec / (dxyz ** 3)
)
if coordinates.lower() == "cylindrical":
b = spatial.cartesian_2_cylindrical(xyz, b)
return b
[docs] def magnetic_field(self, xyz, coordinates="cartesian"):
"""Magnetic field (:math:`\\vec{h}`) of a static magnetic dipole
**Required**
:param numpy.ndarray xyz: Location of the receivers(s)
**Optional**
:param str coordinates: coordinate system that the xyz is provided
in and that the solution will be returned
in (cartesian or cylindrical).
Default: `"cartesian"`
**Returns**
:rtype: numpy.ndarray
:return: The magnetic field at each observation location
"""
return self.magnetic_flux(xyz, coordinates=coordinates) / self.mu
[docs]class CircularLoopWholeSpace(BaseDipole, BaseEM):
"""
Static magnetic field from a circular loop in a wholespace.
"""
current = properties.Float(
"Electric current through the loop (A)", default=1.
)
radius = properties.Float(
"radius of the loop (m)", default=1., min=0.
)
[docs] def vector_potential(self, xyz, coordinates="cartesian"):
"""Vector potential due to the a steady-state current through a
circular loop. We solve in cylindrical coordinates
.. math::
A_\\theta(\\rho, z) = \\frac{\mu_0 I}{\pi k}
\sqrt{R / \\rho^2}[(1 - k^2/2) * K(k^2) - K(k^2)]
where
.. math::
k^2 = \\frac{4 R \\rho}{(R + \\rho)^2 + z^2}
and
- :math:`\\rho = \sqrt{x^2 + y^2}` is the horizontal distance to the test point
- :math:`r` is the distance to a test point
- :math:`I` is the current through the loop
- :math:`R` is the radius of the loop
- :math:`E(k^2)` and :math:`K(k^2)` are the complete elliptic integrals
**Required**
:param numpy.ndarray xyz: Location where we calculate the vector
potential
**Optional**
:param str coordinates: coordinate system that the xyz is provided
in and that the solution will be returned
in (cartesian or cylindrical).
Default: `"cartesian"`
**Returns**
:rtype: numpy.ndarray
:return: The magnetic vector potential at each observation location
"""
eps = 1e-10
supported_coordinates = ["cartesian", "cylindrical"]
assert coordinates.lower() in supported_coordinates, (
"coordinates must be in {}, the coordinate system "
"you provided, {}, is not yet supported".format(
supported_coordinates, coordinates
)
)
# convert coordinates if not cartesian
if coordinates.lower() == "cylindrical":
xyz = spatial.cylindrical_2_cartesian(xyz)
xyz = spatial.rotate_points_from_normals(
xyz, np.array(self.orientation), # work around for a properties issue
np.r_[0., 0., 1.], x0=np.array(self.location)
)
n_obs = xyz.shape[0]
dxyz = self.vector_distance(xyz)
r = self.distance(xyz)
rho = np.sqrt((dxyz[:, :2]**2).sum(1))
k2 = (4 * self.radius * rho) / ((self.radius + rho)**2 +dxyz[:, 2]**2)
k2[k2 > 1.] = 1. # if there are any rounding errors
E = ellipe(k2)
K = ellipk(k2)
# singular if rho = 0, k2 = 1
ind = (rho > eps) & (k2 < 1)
Atheta = np.zeros_like(r)
Atheta[ind] = (
(self.mu * self.current) / (np.pi * np.sqrt(k2[ind])) *
np.sqrt(self.radius / rho[ind]) *
((1. - k2[ind] / 2.)*K[ind] - E[ind])
)
# assume that the z-axis aligns with the polar axis
A = np.zeros_like(xyz)
A[ind, 0] = Atheta[ind] * (-dxyz[ind, 1] / rho[ind])
A[ind, 1] = Atheta[ind] * (dxyz[ind, 0] / rho[ind])
# rotate the points to aligned with the normal to the source
A = spatial.rotate_points_from_normals(
A, np.r_[0., 0., 1.], np.array(self.orientation),
x0=np.array(self.location)
)
if coordinates.lower() == "cylindrical":
A = spatial.cartesian_2_cylindrical(xyz, A)
return A