From bb2c3ec7181ae528e768a2c3a60dc3fd7d8ae44c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 25 Apr 2017 18:39:06 -0700 Subject: [PATCH 01/41] first pass at TDEM dipole --- geoana/em/__init__.py | 1 + geoana/em/base.py | 142 +++++++++++++++++++++++++++++++-------- geoana/em/biot_savart.py | 9 +++ geoana/em/fdem.py | 2 + geoana/em/static.py | 15 +++-- geoana/em/tdem.py | 89 ++++++++++++++++++++++++ geoana/spatial.py | 36 ++++++++++ 7 files changed, 260 insertions(+), 34 deletions(-) create mode 100644 geoana/em/biot_savart.py create mode 100644 geoana/em/tdem.py create mode 100644 geoana/spatial.py diff --git a/geoana/em/__init__.py b/geoana/em/__init__.py index fda393e5..f733b529 100644 --- a/geoana/em/__init__.py +++ b/geoana/em/__init__.py @@ -1,2 +1,3 @@ from . import static from . import fdem +from . import tdem diff --git a/geoana/em/base.py b/geoana/em/base.py index 9dcc58f4..9f2cc00a 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -7,23 +7,100 @@ import properties from scipy.constants import mu_0, pi, epsilon_0 +from .. import spatial + + +############################################################################### +# # +# Functions # +# # +############################################################################### + +def omega(frequency): + """ + Angular frequency + + :param frequency float: frequency (Hz) + """ + return 2*np.pi*frequency + + +def wave_number(frequency, sigma, mu=mu_0, epsilon=epsilon_0): + """ + Wavenumber of an electromagnetic wave in a medium with constant physical + properties + + :param frequency float: frequency (Hz) + :param sigma float: electrical conductivity (S/m) + :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m + :param epsilon float: dielectric permittivity (F/m). Default: :math:`\epsilon_0 = 8.85 \times 10^{-12}` F/m + """ + omega = omega(frequency) + return np.sqrt(omega**2. * mu * epsilon - 1j * omega * mu * sigma) + + +def skin_depth(frequency, sigma, mu=mu_0): + """ + Distance at which an em wave has decayed by a factor of 1/e in a medium + with constant physical properties + + :param frequency float: frequency (Hz) + :param sigma float: electrical conductivity (S/m) + :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m + """ + omega = omega(frequency) + return np.sqrt(2./(omega*sigma*mu)) + + +def peak_time(z, sigma, mu=mu_0): + """ + Time at which the maximum signal amplitude is observed at a particular + location for a transient plane wave through a homogeneous medium. + + See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html + + :param z float: distance from source (m) + :param sigma float: electrical conductivity (S/m) + :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m + """ + return (mu * sigma * z**2)/6. + + +def diffusion_distance(time, sigma, mu=mu_0): + """ + Distance at which the signal amplitude is largest for a given time after + shut off. Also referred to as the peak distance + + See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html + + + """ + return np.sqrt(2*time/(mu*sigma)) + + +############################################################################### +# # +# Base Classes # +# # +############################################################################### + class BaseEM(properties.HasProperties): mu = properties.Float( - help='Magnetic permeability.', + "Magnetic permeability.", default=mu_0, min=0.0 ) sigma = properties.Float( - help='Electrical conductivity (S/m)', + "Electrical conductivity (S/m)", default=1.0, min=0.0 ) epsilon = properties.Float( - help='Permitivity value', + "Permitivity value", default=epsilon_0, min=0.0 ) @@ -32,42 +109,34 @@ class BaseEM(properties.HasProperties): class BaseDipole(BaseEM): orientation = properties.Vector3( - help='orientation of dipole', - default='X', + "orientation of dipole", + default="X", length=1.0 ) location = properties.Vector3( - help='location of the electric dipole source', - default='ZERO' + "location of the electric dipole source", + default="ZERO" ) - def offset_from_location(self, xyz): + def vector_distance(self, xyz): + return spatial.vector_distance(xyz, self.location) - # TODO: validate stuff - # xyz = Utils.asArray_N_x_Dim(xyz, 3) - - return np.c_[ - xyz[:, 0] - self.location[0], - xyz[:, 1] - self.location[1], - xyz[:, 2] - self.location[2] - ] - - def distance_from_location(self, xyz): - return np.sqrt((self.offset_from_location(xyz)**2).sum(axis=1)) + def distance(self, xyz): + return spatial.distance(xyz, self.location) class BaseFDEM(BaseEM): frequency = properties.Float( - help='Source frequency (Hz)', + "Source frequency (Hz)", default=1e2, min=0.0 ) @property def omega(self): - return 2.0*pi*self.frequency + return omega(self.frequency) @property def sigma_hat(self): @@ -75,22 +144,39 @@ def sigma_hat(self): @property def wave_number(self): - np.sqrt( - self.omega**2. * self.mu * self.epsilon - - 1j * self.omega * self.mu * self.sigma - ) + return wave_number(self.frequency, self.sigma, self.mu) + + @property + def skin_depth(self): + return skin_depth(self.frequency, self.sigma, self.mu) + + +class BaseTDEM(BaseEM): + + time = properties.Float( + "time after shut-off at which we are evaluating the fields (s)", + + required=True + ) + + 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) class BaseElectricDipole(BaseDipole): length = properties.Float( - help='length of the dipole (m)', + "length of the dipole (m)", default=1.0, min=0.0 ) current = properties.Float( - help='size of the injected current (A)', + "size of the injected current (A)", default=1.0, min=0.0 ) @@ -99,7 +185,7 @@ class BaseElectricDipole(BaseDipole): class BaseMagneticDipole(BaseDipole): moment = properties.Float( - help='moment of the dipole (Am^2)', + "moment of the dipole (Am^2)", default=1.0, min=0.0 ) diff --git a/geoana/em/biot_savart.py b/geoana/em/biot_savart.py new file mode 100644 index 00000000..e4c8c062 --- /dev/null +++ b/geoana/em/biot_savart.py @@ -0,0 +1,9 @@ +import numpy as np +import scipy.sparse as sp + +from SimPEG import Mesh, Utils + +from .base import BaseEM + +class BiotSavart(BaseEM): + pass diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index aef7784b..24aef7c9 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -16,6 +16,8 @@ def electric_field(self, xyz, **kwargs): pass + + def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): """E_from_EDWS Computing the analytic electric fields (E) from an electrical dipole in a wholespace diff --git a/geoana/em/static.py b/geoana/em/static.py index 2f76e38b..64fecd65 100644 --- a/geoana/em/static.py +++ b/geoana/em/static.py @@ -28,9 +28,9 @@ def vector_potential(self, xyz, **kwargs): offset = self.offset_from_location(xyz) dist = self.distance_from_location(xyz) - m_vec = self.moment * np.atleast_2d( - self.orientation - ).repeat(n_obs, axis=0) + m_vec = ( + self.moment * np.atleast_2d(self.orientation).repeat(n_obs, axis=0) + ) # Repeat the scalars dist = np.atleast_2d(dist).T.repeat(3, axis=1) @@ -52,9 +52,9 @@ def magnetic_flux(self, xyz, **kwargs): offset = self.offset_from_location(xyz) dist = self.distance_from_location(xyz) - m_vec = self.moment * np.atleast_2d( - self.orientation - ).repeat(n_obs, axis=0) + m_vec = ( + self.moment * np.atleast_2d(self.orientation).repeat(n_obs, axis=0) + ) m_dot_r = (m_vec * offset).sum(axis=1) @@ -77,3 +77,6 @@ def magnetic_flux_equation(): def magnetic_field(self, xyz, **kwargs): return self.magnetic_flux(xyz, **kwargs) / self.mu + + +# class diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py new file mode 100644 index 00000000..1730a622 --- /dev/null +++ b/geoana/em/tdem.py @@ -0,0 +1,89 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + +from .base import BaseElectricDipole, BaseTDEM + +from scipy.constants import mu_0, pi, epsilon_0 +from scipy.special import erf +import numpy as np +import warnings + + +class ElectricDipole(BaseTDEM, BaseElectricDipole): + + """ + .. todo: write this independent of source orientation with dot products + """ + @property + def theta(self): + return np.sqrt(self.mu*self.sigma / (4.*self.time)) + + def electric_field(self, xyz): + dxyz = self.vector_distance(xyz) + x = dxyz[:, 0] + y = dxyz[:, 1] + z = dxyz[:, 2] + + root_pi = np.sqrt(np.pi) + r = self.distance(xyz) + r2 = r**2 + r3 = r**3. + + theta_r = self.theta * r + e_n_theta_r_2 = np.exp(-theta_r**2) + + erf_thetat_r = erf(theta_r) + + current_term = ( + (self.current * self.length) / + (4. * np.pi * self.sigma * r3) + ) + + symmetric_term = ( + - ( + (4./root_pi) * theta_r**3 + (6./root_pi) * theta_r + ) * e_n_theta_r_2 + + 3 * erf_thetat_r + ) + + src_orientation_term = ( + - ( + (4./root_pi) * theta_r**3 + (2./root_pi) * theta_r + ) * e_n_theta_r_2 + + erf_thetat_r + ) + + if np.all(self.orientation == np.r_[1., 0., 0.]): + ex = current_term * ( + symmetric_term * (x**2/r2) - src_orientation_term + ) + ey = current_term * ( symmetric_term * (x*y)/r2 ) + ez = current_term * ( symmetric_term * (x*z)/r2 ) + + elif np.all(self.orientation == np.r_[0., 1., 0.]): + ey = current_term * ( + symmetric_term * (y**2/r2) - src_orientation_term + ) + ez = current_term * ( symmetric_term * (y*z)/r2 ) + ex = current_term * ( symmetric_term * (y*x)/r2 ) + + elif np.all(self.orientation == np.r_[0., 0., 1.]): + ez = current_term * ( + symmetric_term * (z**2/r2) - src_orientation_term + ) + ex = current_term * ( symmetric_term * (z*x)/r2 ) + ey = current_term * ( symmetric_term * (z*y)/r2 ) + else: + raise NotImplementedError + + return np.c_[ex, ey, ez] + + def current_density(self, xyz): + return self.sigma * self.e(xyz) + + + + + diff --git a/geoana/spatial.py b/geoana/spatial.py new file mode 100644 index 00000000..559a3107 --- /dev/null +++ b/geoana/spatial.py @@ -0,0 +1,36 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + +import numpy as np +import properties + + +def vector_distance(xyz, origin=np.r_[0., 0., 0.]): + """ + Vector distance of a grid, xyz from an origin origin. + :param numpy.ndarray xyz: grid (npoints x 3) + :param numpy.ndarray origin: origin (default: [0., 0., 0.]) + """ + assert(xyz.shape[1] == 3), ( + "the xyz grid should be npoints by 3, the shape provided is {}".format( + xyz.shape + ) + ) + + dx = xyz[:, 0] - origin[0] + dy = xyz[:, 1] - origin[1] + dz = xyz[:, 2] - origin[2] + + return np.c_[dx, dy, dz] + + +def distance(xyz, origin=np.r_[0., 0., 0.]): + """ + Radial distance from an origin origin + :param numpy.ndarray xyz: grid (npoints x 3) + :param numpy.ndarray origin: origin (default: [0., 0., 0.]) + """ + dxyz = vector_distance(xyz, origin) + return np.sqrt((dxyz**2).sum(axis=1)) From 776f743f4067d2a726b1f7939e6f4252e8d0dc55 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 6 May 2017 14:20:22 -0700 Subject: [PATCH 02/41] minor cleanup in TDEM dipoles --- geoana/em/tdem.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 1730a622..8b0e70ca 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -16,6 +16,9 @@ class ElectricDipole(BaseTDEM, BaseElectricDipole): """ .. todo: write this independent of source orientation with dot products """ + def __init__(self, **kwargs): + super(ElectricDipole, self).__init__(**kwargs) + @property def theta(self): return np.sqrt(self.mu*self.sigma / (4.*self.time)) @@ -81,9 +84,5 @@ def electric_field(self, xyz): return np.c_[ex, ey, ez] def current_density(self, xyz): - return self.sigma * self.e(xyz) - - - - + return self.sigma * self.electric_field(xyz) From dafbbbdef68f4d5dcb1fe926432abe44b794d944 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 22 Jun 2017 16:42:52 -0700 Subject: [PATCH 03/41] add docs for em base classes, move baseFDEM to fdem. --- geoana/em/base.py | 53 ++++++++++++++++++--------------------------- geoana/em/fdem.py | 51 ++++++++++++++++++++++++++++++++++++++++--- geoana/em/static.py | 4 ++-- 3 files changed, 71 insertions(+), 37 deletions(-) diff --git a/geoana/em/base.py b/geoana/em/base.py index 9dcc58f4..16f33729 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -9,36 +9,43 @@ class BaseEM(properties.HasProperties): + """ + Base class for electromanetics. Contains physical properties that are + relevant to all problems that use Maxwell's equations + """ mu = properties.Float( - help='Magnetic permeability.', + "Magnetic permeability (H/m)", default=mu_0, min=0.0 ) sigma = properties.Float( - help='Electrical conductivity (S/m)', + "Electrical conductivity (S/m)", default=1.0, min=0.0 ) epsilon = properties.Float( - help='Permitivity value', + "Permitivity value (F/m)", default=epsilon_0, min=0.0 ) class BaseDipole(BaseEM): + """ + Base class for dipoles. + """ orientation = properties.Vector3( - help='orientation of dipole', + "orientation of dipole", default='X', length=1.0 ) location = properties.Vector3( - help='location of the electric dipole source', + "location of the electric dipole source", default='ZERO' ) @@ -57,49 +64,31 @@ def distance_from_location(self, xyz): return np.sqrt((self.offset_from_location(xyz)**2).sum(axis=1)) -class BaseFDEM(BaseEM): - - frequency = properties.Float( - help='Source frequency (Hz)', - default=1e2, - min=0.0 - ) - - @property - def omega(self): - return 2.0*pi*self.frequency - - @property - def sigma_hat(self): - return self.sigma + 1j*self.omega*self.epsilon - - @property - def wave_number(self): - np.sqrt( - self.omega**2. * self.mu * self.epsilon - - 1j * self.omega * self.mu * self.sigma - ) - - class BaseElectricDipole(BaseDipole): + """ + Base class for electric current dipoles + """ length = properties.Float( - help='length of the dipole (m)', + "length of the dipole (m)", default=1.0, min=0.0 ) current = properties.Float( - help='size of the injected current (A)', + "size of the injected current (A)", default=1.0, min=0.0 ) class BaseMagneticDipole(BaseDipole): + """ + Base class for magnetic dipoles + """ moment = properties.Float( - help='moment of the dipole (Am^2)', + "moment of the dipole (Am^2)", default=1.0, min=0.0 ) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index aef7784b..f4e3ad11 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -3,18 +3,63 @@ from __future__ import print_function from __future__ import unicode_literals -from .base import BaseElectricDipole, BaseFDEM - from scipy.constants import mu_0, pi, epsilon_0 import numpy as np import warnings +from .base import BaseElectricDipole, BaseMagneticDipole + + +class BaseFDEM(BaseEM): + """ + Base frequency domain EM class. Contains methds and properties usefull + across all FDEM problems. + """ + + frequency = properties.Float( + "Source frequency (Hz)", + default=1e2, + min=0.0 + ) + + @property + def omega(self): + """ + """ + return 2.0*pi*self.frequency + + @property + def sigma_hat(self): + return self.sigma + 1j*self.omega*self.epsilon + + @property + def wave_number(self): + return np.sqrt( + self.omega**2. * self.mu * self.epsilon - + 1j * self.omega * self.mu * self.sigma + ) + + @property + def skin_depth(self): + return np.sqrt() + -class ElectricDipole_WholeSpace(BaseElectricDipole, BaseFDEM): +class ElectricDipoleWholeSpace(BaseElectricDipole, BaseFDEM): def electric_field(self, xyz, **kwargs): pass + def current_density(self, xyz, **kwargs): + pass + + def magnetic_field(self, xyz, **kwargs): + pass + + def magnetic_flux_density(self, xyz, **kwargs): + pass + + +class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM) def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): """E_from_EDWS diff --git a/geoana/em/static.py b/geoana/em/static.py index 2f76e38b..5db0cda2 100644 --- a/geoana/em/static.py +++ b/geoana/em/static.py @@ -3,12 +3,12 @@ from __future__ import print_function from __future__ import unicode_literals -from .base import BaseEM, BaseMagneticDipole - import numpy as np from ipywidgets import Latex import warnings +from .base import BaseEM, BaseMagneticDipole, BaseElectricDipole + class MagneticDipole_WholeSpace(BaseMagneticDipole, BaseEM): From 512079898374fd43f35a61e96d48968177a94000 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 22 Jun 2017 16:45:53 -0700 Subject: [PATCH 04/41] remove biot savart --- geoana/em/biot_savart.py | 9 --------- 1 file changed, 9 deletions(-) delete mode 100644 geoana/em/biot_savart.py diff --git a/geoana/em/biot_savart.py b/geoana/em/biot_savart.py deleted file mode 100644 index e4c8c062..00000000 --- a/geoana/em/biot_savart.py +++ /dev/null @@ -1,9 +0,0 @@ -import numpy as np -import scipy.sparse as sp - -from SimPEG import Mesh, Utils - -from .base import BaseEM - -class BiotSavart(BaseEM): - pass From f5a1d350299f696fa02cf7871b505b9ee7ccdff6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 22 Jun 2017 17:57:33 -0700 Subject: [PATCH 05/41] documentation and organizing --- geoana/em/base.py | 66 -------------------- geoana/em/fdem.py | 156 ++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 145 insertions(+), 77 deletions(-) diff --git a/geoana/em/base.py b/geoana/em/base.py index 88d9f1a2..671a2ffc 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -10,48 +10,6 @@ from .. import spatial -############################################################################### -# # -# Functions # -# # -############################################################################### - -def omega(frequency): - """ - Angular frequency - - :param frequency float: frequency (Hz) - """ - return 2*np.pi*frequency - - -def wave_number(frequency, sigma, mu=mu_0, epsilon=epsilon_0): - """ - Wavenumber of an electromagnetic wave in a medium with constant physical - properties - - :param frequency float: frequency (Hz) - :param sigma float: electrical conductivity (S/m) - :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m - :param epsilon float: dielectric permittivity (F/m). Default: :math:`\epsilon_0 = 8.85 \times 10^{-12}` F/m - """ - omega = omega(frequency) - return np.sqrt(omega**2. * mu * epsilon - 1j * omega * mu * sigma) - - -def skin_depth(frequency, sigma, mu=mu_0): - """ - Distance at which an em wave has decayed by a factor of 1/e in a medium - with constant physical properties - - :param frequency float: frequency (Hz) - :param sigma float: electrical conductivity (S/m) - :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m - """ - omega = omega(frequency) - return np.sqrt(2./(omega*sigma*mu)) - - def peak_time(z, sigma, mu=mu_0): """ Time at which the maximum signal amplitude is observed at a particular @@ -133,30 +91,6 @@ def distance(self, xyz): return spatial.distance(xyz, self.location) -class BaseFDEM(BaseEM): - - frequency = properties.Float( - "Source frequency (Hz)", - default=1e2, - min=0.0 - ) - - @property - def omega(self): - return omega(self.frequency) - - @property - def sigma_hat(self): - return self.sigma + 1j*self.omega*self.epsilon - - @property - def wave_number(self): - return wave_number(self.frequency, self.sigma, self.mu) - - @property - def skin_depth(self): - return skin_depth(self.frequency, self.sigma, self.mu) - class BaseTDEM(BaseEM): diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 3fc88cb2..06d964e1 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -7,44 +7,167 @@ import numpy as np import warnings -from .base import BaseElectricDipole, BaseMagneticDipole +from .base import BaseElectricDipole, BaseMagneticDipole, BaseEM -class BaseFDEM(BaseEM): +############################################################################### +# # +# Utility Functions # +# # +############################################################################### + +def omega(frequency): + """ + Angular frequency + + .. math:: + + \omega = 2 \pi f + + :param frequency float: frequency (Hz) + """ + return 2*np.pi*frequency + + +def wave_number( + frequency, sigma, mu=mu_0, epsilon=epsilon_0, quasi_static=False +): + """ + Wavenumber of an electromagnetic wave in a medium with constant physical + properties + + .. math:: + + k = \sqrt{\omega**2 \mu \varepsilon - i \omega \mu \sigma} + + :param (float, numpy.ndarray) frequency: frequency (Hz) + :param float sigma: electrical conductivity (S/m) + :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 quasi_static: use the quasi-static assumption? Default: False + """ + omega = omega(frequency) + if quasi_static is True: + return np.sqrt(-1j * omega * mu * sigma) + return np.sqrt(omega**2. * mu * epsilon - 1j * omega * mu * sigma) + + +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}} + + :param float frequency: frequency (Hz) + :param float sigma: electrical conductivity (S/m) + :param float mu: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m """ - Base frequency domain EM class. Contains methds and properties usefull - across all FDEM problems. + omega = omega(frequency) + return np.sqrt(2./(omega*sigma*mu)) + + +def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasi_static=False): + """ + conductivity with displacement current contribution + + .. math:: + + \hat{\sigma} = \sigma + i \omega \varepsilon + + :param (float, numpy.array) frequency: frequency (Hz) + :param float sigma: electrical conductivity (S/m) + :param float epsilon: dielectric permittivity. Default :math:`\varepsilon_0` + :param bool quasi_static: use the quasi-static assumption? Default: False """ + if quasi_static is True: + return sigma + return sigma + 1j*omega(frequency)*epsilon + +############################################################################### +# # +# Classes # +# # +############################################################################### + +class BaseFDEM(BaseEM): + """ + Base frequency domain electromagnetic class + """ frequency = properties.Float( "Source frequency (Hz)", default=1e2, min=0.0 ) + quasi_static = 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 2.0*pi*self.frequency + return omega(self.frequency) @property def sigma_hat(self): - return self.sigma + 1j*self.omega*self.epsilon + """ + conductivity with displacement current contribution + + .. math:: + + \hat{\sigma} = \sigma + i \omega \varepsilon + + """ + return sigma_hat( + self.frequency, self.sigma, self.epsilon, self.quasi_static + ) @property def wave_number(self): - return np.sqrt( - self.omega**2. * self.mu * self.epsilon - - 1j * self.omega * self.mu * self.sigma - ) + """ + Wavenumber of an electromagnetic wave in a medium with constant physical + properties + + .. math:: + + k = \sqrt{\omega**2 \mu \varepsilon - i \omega \mu \sigma} """ + return wave_number(self.frequency, self.sigma, self.mu) @property def skin_depth(self): - return np.sqrt() + """ + 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, self.mu) 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:`\hat{\mathbf{r}}_s` + + .. math:: + + \mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r} - \mathbf{r}_s) + """ def electric_field(self, xyz, **kwargs): pass @@ -61,6 +184,17 @@ def magnetic_flux_density(self, xyz, **kwargs): class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM) + def electric_field(self, xyz, **kwargs): + pass + + def current_density(self, xyz, **kwargs): + pass + + def magnetic_field(self, xyz, **kwargs): + pass + + def magnetic_flux_density(self, xyz, **kwargs): + pass def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): From 785113719bb51998f2785eeae076e1c564521163 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 23 Jun 2017 00:18:58 -0700 Subject: [PATCH 06/41] electric dipole in a wholespace initial implementation --- geoana/em/base.py | 21 ++++++++++ geoana/em/fdem.py | 102 ++++++++++++++++++++++++++++++++++++++-------- geoana/spatial.py | 20 +++++++-- 3 files changed, 123 insertions(+), 20 deletions(-) diff --git a/geoana/em/base.py b/geoana/em/base.py index 671a2ffc..231a2872 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -85,11 +85,32 @@ class BaseDipole(BaseEM): ) def vector_distance(self, xyz): + """ + Vector distance from the dipole location + :param numpy.ndarray xyz: grid + """ return spatial.vector_distance(xyz, self.location) def distance(self, xyz): + """ + Distance from the dipole location + """ return spatial.distance(xyz, self.location) + def dot_orientation(self, xyz): + """ + Take the dot product between a grid and the orientation of the dipole + """ + return spatial.vector_dot(xyz, self.orientation) + + def cross_orientation(self, xyz): + """ + Take the cross product between a grid and the orientation of the dipole + """ + orientation = np.kron( + np.atleast_2d(self.orientation), np.ones(xyz.shape[0]).T + ) + return np.cross(xyz, orientation) class BaseTDEM(BaseEM): diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 06d964e1..e6ce4a9e 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -157,43 +157,113 @@ def skin_depth(self): return skin_depth(self.frequency, self.sigma, self.mu) -class ElectricDipoleWholeSpace(BaseElectricDipole, BaseFDEM): +class BaseFDEMDipoleWholeSpace(BaseFDEM): + """ + Base FDEM Dipole + """ + + +class ElectricDipoleWholeSpace( + BaseElectricDipole, BaseFDEMDipoleWholeSpace +): """ 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:`\hat{\mathbf{r}}_s` + :math:`\mathbf{\hat{u}}` .. math:: - \mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r} - \mathbf{r}_s) + \mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r} - \mathbf{r}_s)\mathbf{\hat{u}} """ + def vector_potential(self, xyz): + """ + Vector potential for an electric dipole in a wholespace - def electric_field(self, xyz, **kwargs): - pass + .. math:: - def current_density(self, xyz, **kwargs): - pass + \mathbf{A} = \frac{I ds}{4 \pi r} e^{-ikr}\mathbf{\hat{u}} - def magnetic_field(self, xyz, **kwargs): - pass + """ + r = self.distance(xyz) + a = ( + (self.current * self.length) / (4*np.pi*r) * + np.exp(-i*self.wave_number*r) + ) + a = np.kron(np.ones(1, 3), np.atleast_2d(a).T) + return self.dot_orientation(a) - def magnetic_flux_density(self, xyz, **kwargs): - pass + def electric_field(self, xyz): + """ + Electric field from an electric dipole + .. math:: -class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM) + \mathbf{E} = \frac{1}{\hat{\sigma}} \nabla \nabla \cdot \mathbf{A} - i \omega \mu \mathbf{A} - def electric_field(self, xyz, **kwargs): + """ + dxyz = self.vector_distance(xyz) + r = self.distance(xyz) + r = np.kron(np.ones(1, 3), np.atleast_2d(r).T) + kr = self.wave_number * r + + front = ( + (self.current * self.length) / (4 * np.pi * self.sigma * r**3 ) * + np.exp(-1j * kr) + ) + symmetric_term = ( + self.dot_orientation(dxyz) * (-kr**2 + 3*1j*kr + 3) / r**2 + ) + oriented_term = ( + (kr**2 - 1j*kr - 1) * + np.kron(self.orientation, np.ones(dxyz.shape[0], 1)) + ) + return front * ( symmetric_term + oriented_term ) + + def current_density(self, xyz): + return self.sigma * self.electric_field(xyz) + + def magnetic_field(self, xyz): + """ + Magnetic field from an electric dipole + + .. math:: + + \mathbf{H} = \nabla \times \mathbf{A} + + """ + dxyz = self.vector_distance(xyz) + r = self.distance(xyz) + r = np.kron(np.ones(1, 3), np.atleast_2d(r).T) + kr = self.wave_number * r + + front = ( + self.current * self.length / (4 * np.pi * r**2) * (1j * kr + 1) * + np.exp(-1j * kr) + ) + return front * self.cross_orientation(xyz) / r + + def magnetic_flux_density(self, xyz): + """ + magnetic flux density from an electric dipole + """ + return self.mu * self.magnetic_field(xyz) + + +class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM) + """ + Harmonic magnetic dipole in a whole space. + """ + def electric_field(self, xyz): pass - def current_density(self, xyz, **kwargs): + def current_density(self, xyz): pass - def magnetic_field(self, xyz, **kwargs): + def magnetic_field(self, xyz): pass - def magnetic_flux_density(self, xyz, **kwargs): + def magnetic_flux_density(self, xyz): pass diff --git a/geoana/spatial.py b/geoana/spatial.py index 559a3107..d6d5f99c 100644 --- a/geoana/spatial.py +++ b/geoana/spatial.py @@ -10,8 +10,8 @@ def vector_distance(xyz, origin=np.r_[0., 0., 0.]): """ Vector distance of a grid, xyz from an origin origin. - :param numpy.ndarray xyz: grid (npoints x 3) - :param numpy.ndarray origin: origin (default: [0., 0., 0.]) + :param numpy.array xyz: grid (npoints x 3) + :param numpy.array origin: origin (default: [0., 0., 0.]) """ assert(xyz.shape[1] == 3), ( "the xyz grid should be npoints by 3, the shape provided is {}".format( @@ -29,8 +29,20 @@ def vector_distance(xyz, origin=np.r_[0., 0., 0.]): def distance(xyz, origin=np.r_[0., 0., 0.]): """ Radial distance from an origin origin - :param numpy.ndarray xyz: grid (npoints x 3) - :param numpy.ndarray origin: origin (default: [0., 0., 0.]) + :param numpy.array xyz: grid (npoints x 3) + :param numpy.array origin: origin (default: [0., 0., 0.]) """ dxyz = vector_distance(xyz, origin) return np.sqrt((dxyz**2).sum(axis=1)) + + +def vector_dot(xyz, vector): + """ + Take a dot product between an array of vectors, xyz and a vector [x, y, z] + :param numpy.array xyz: grid (npoints x 3) + :param numpy.array vector: vector (1 x 3) + """ + assert len(vector) == 3, "vector should be length 3" + return np.hstack([ + vector[0]*xyz[:, 0], vector[1]*xyz[:, 1], vector[2]*xyz[:, 2] + ]) From dedf6a5d16f3d386b513c54029cd211b1211edd6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 25 Jun 2017 16:06:06 -0700 Subject: [PATCH 07/41] bug fixes in imports, start docs --- docs/content/earthquake.rst | 8 ++++++ docs/content/em.rst | 8 ++++++ docs/index.rst | 2 ++ geoana/em/fdem.py | 56 ++----------------------------------- geoana/em/static.py | 3 +- tests/test_em_fdem.py | 56 +++++++++++++++++++++++++++++++++++++ 6 files changed, 78 insertions(+), 55 deletions(-) create mode 100644 docs/content/earthquake.rst create mode 100644 docs/content/em.rst diff --git a/docs/content/earthquake.rst b/docs/content/earthquake.rst new file mode 100644 index 00000000..396dedbd --- /dev/null +++ b/docs/content/earthquake.rst @@ -0,0 +1,8 @@ +.. _earthquake: + +Earthquake +========== + +.. autoclass:: geoana.earthquake + :members: + :undoc-members: diff --git a/docs/content/em.rst b/docs/content/em.rst new file mode 100644 index 00000000..5c65476d --- /dev/null +++ b/docs/content/em.rst @@ -0,0 +1,8 @@ +.. _em: + +Electromagnetics +================ + +.. autoclass:: geoana.em + :members: + :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index 2c7c7116..50ec6668 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,6 +7,8 @@ Contents: .. toctree:: :maxdepth: 2 + content/earthquake + content/em Indices and tables diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index e6ce4a9e..347b97b2 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -6,6 +6,7 @@ from scipy.constants import mu_0, pi, epsilon_0 import numpy as np import warnings +import properties from .base import BaseElectricDipole, BaseMagneticDipole, BaseEM @@ -161,6 +162,7 @@ class BaseFDEMDipoleWholeSpace(BaseFDEM): """ Base FDEM Dipole """ + pass class ElectricDipoleWholeSpace( @@ -250,7 +252,7 @@ def magnetic_flux_density(self, xyz): return self.mu * self.magnetic_field(xyz) -class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM) +class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM): """ Harmonic magnetic dipole in a whole space. """ @@ -267,58 +269,6 @@ def magnetic_flux_density(self, xyz): pass -def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): - """E_from_EDWS - Computing the analytic electric fields (E) from an electrical dipole in a wholespace - - You have the option of computing E for multiple frequencies at a single reciever location - or a single frequency at multiple locations - - :param numpy.array XYZ: reciever locations at which to evaluate E - :param float epsr: relative permitivitty value (unitless), default is 1.0 - :rtype: numpy.array - :return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and frequencies. - """ - - mu = mu_0*(1+kappa) - epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omega(f)*epsilon - - XYZ = Utils.asArray_N_x_Dim(XYZ, 3) - # Check - if XYZ.shape[0] > 1 & f.shape[0] > 1: - raise Exception('I/O type error: For multiple field locations only a single frequency can be specified.') - - dx = XYZ[:, 0] - srcLoc[0] - dy = XYZ[:, 1] - srcLoc[1] - dz = XYZ[:, 2] - srcLoc[2] - - r = np.sqrt(dx**2. + dy**2. + dz**2.) - # k = np.sqrt( -1j*2.*pi*f*mu*sig ) - k = np.sqrt(omega(f)**2. * mu * epsilon - 1j * omega(f) * mu * sig) - - front = current * length / (4.*pi*sig_hat* r**3) * np.exp(-1j*k*r) - mid = -k**2 * r**2 + 3*1j*k*r + 3 - - if orientation.upper() == 'X': - Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ey = front*(dx*dy / r**2)*mid - Ez = front*(dx*dz / r**2)*mid - return Ex, Ey, Ez - - elif orientation.upper() == 'Y': - # x--> y, y--> z, z-->x - Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ez = front*(dy*dz / r**2)*mid - Ex = front*(dy*dx / r**2)*mid - return Ex, Ey, Ez - - elif orientation.upper() == 'Z': - # x --> z, y --> x, z --> y - Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ex = front*(dz*dx / r**2)*mid - Ey = front*(dz*dy / r**2)*mid - return Ex, Ey, Ez - def MagneticDipoleFields(srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0): """ diff --git a/geoana/em/static.py b/geoana/em/static.py index 2106564e..125f26d1 100644 --- a/geoana/em/static.py +++ b/geoana/em/static.py @@ -4,7 +4,6 @@ from __future__ import unicode_literals import numpy as np -from ipywidgets import Latex import warnings from .base import BaseEM, BaseMagneticDipole, BaseElectricDipole @@ -70,7 +69,7 @@ def magnetic_flux(self, xyz, **kwargs): @staticmethod def magnetic_flux_equation(): - return Latex( + return "{}".format( "$\\frac{\mu}{4\pi} \\frac{\mathbf{m} \\times " "\mathbf{\hat{r}}}{r^2}$" ) diff --git a/tests/test_em_fdem.py b/tests/test_em_fdem.py index 876df2c1..987fc862 100644 --- a/tests/test_em_fdem.py +++ b/tests/test_em_fdem.py @@ -8,10 +8,66 @@ from geoana.em import fdem +def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): + """E_from_EDWS + Computing the analytic electric fields (E) from an electrical dipole in a wholespace + - You have the option of computing E for multiple frequencies at a single reciever location + or a single frequency at multiple locations + + :param numpy.array XYZ: reciever locations at which to evaluate E + :param float epsr: relative permitivitty value (unitless), default is 1.0 + :rtype: numpy.array + :return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and frequencies. + """ + + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omega(f)*epsilon + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception('I/O type error: For multiple field locations only a single frequency can be specified.') + + dx = XYZ[:, 0] - srcLoc[0] + dy = XYZ[:, 1] - srcLoc[1] + dz = XYZ[:, 2] - srcLoc[2] + + r = np.sqrt(dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*pi*f*mu*sig ) + k = np.sqrt(omega(f)**2. * mu * epsilon - 1j * omega(f) * mu * sig) + + front = current * length / (4.*pi*sig_hat* r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ez = front*(dy*dz / r**2)*mid + Ex = front*(dy*dx / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + + class TestFDEM(unittest.TestCase): def test_vector(self): edws = fdem.ElectricDipole_WholeSpace() + + + if __name__ == '__main__': unittest.main() From 04707dcc39d4f9f602b74001517017df91be697f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 26 Jun 2017 13:52:59 -0700 Subject: [PATCH 08/41] arbitrarily oriented electric dipole --- .travis.yml | 2 +- geoana/em/fdem.py | 118 ++-------------- geoana/spatial.py | 18 ++- requirements_dev.txt | 2 + tests/test_em_fdem.py | 73 ---------- tests/test_em_fdem_electric_dipole.py | 187 ++++++++++++++++++++++++++ tests/test_em_fdem_magnetic_dipole.py | 112 +++++++++++++++ 7 files changed, 326 insertions(+), 186 deletions(-) create mode 100644 requirements_dev.txt delete mode 100644 tests/test_em_fdem.py create mode 100644 tests/test_em_fdem_electric_dipole.py create mode 100644 tests/test_em_fdem_magnetic_dipole.py diff --git a/.travis.yml b/.travis.yml index f610782f..a381b053 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,7 +11,7 @@ env: install: - pip install --upgrade pip; - pip install -r requirements.txt - - pip install sphinx_rtd_theme + - pip install -r requirements_dev.txt - git clone https://github.com/3ptscience/vectormath; - export PYTHONPATH=$PWD/geoana:$PWD/vectormath:$PYTHONPATH diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 347b97b2..7df35295 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -9,6 +9,7 @@ import properties from .base import BaseElectricDipole, BaseMagneticDipole, BaseEM +from .. import spatial ############################################################################### @@ -47,10 +48,10 @@ def wave_number( :param float epsilon: dielectric permittivity (F/m). Default: :math:`\epsilon_0 = 8.85 \times 10^{-12}` F/m :param bool quasi_static: use the quasi-static assumption? Default: False """ - omega = omega(frequency) + w = omega(frequency) if quasi_static is True: - return np.sqrt(-1j * omega * mu * sigma) - return np.sqrt(omega**2. * mu * epsilon - 1j * omega * mu * sigma) + return np.sqrt(-1j * w * mu * sigma) + return np.sqrt(w**2. * mu * epsilon - 1j * w * mu * sigma) def skin_depth(frequency, sigma, mu=mu_0): @@ -100,7 +101,7 @@ class BaseFDEM(BaseEM): """ frequency = properties.Float( "Source frequency (Hz)", - default=1e2, + default=1, min=0.0 ) @@ -206,7 +207,7 @@ def electric_field(self, xyz): """ dxyz = self.vector_distance(xyz) r = self.distance(xyz) - r = np.kron(np.ones(1, 3), np.atleast_2d(r).T) + r = spatial.repeat_scalar(r) kr = self.wave_number * r front = ( @@ -214,11 +215,12 @@ def electric_field(self, xyz): np.exp(-1j * kr) ) symmetric_term = ( - self.dot_orientation(dxyz) * (-kr**2 + 3*1j*kr + 3) / r**2 + spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz * + (-kr**2 + 3*1j*kr + 3) / r**2 ) oriented_term = ( (kr**2 - 1j*kr - 1) * - np.kron(self.orientation, np.ones(dxyz.shape[0], 1)) + np.kron(self.orientation, np.ones((dxyz.shape[0], 1))) ) return front * ( symmetric_term + oriented_term ) @@ -267,105 +269,3 @@ def magnetic_field(self, xyz): def magnetic_flux_density(self, xyz): pass - - - -def MagneticDipoleFields(srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0): - """ - Calculate the vector potential of a set of magnetic dipoles - at given locations 'ref. ' - - .. math:: - - B = \frac{\mu_0}{4 \pi r^3} \left( \frac{3 \vec{r} (\vec{m} \cdot - \vec{r})}{r^2}) - - \vec{m} - \right) \cdot{\hat{rx}} - - :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) - :param numpy.ndarray obsLoc: Where the potentials will be calculated - (x, y, z) - :param str component: The component to calculate - 'x', 'y', or 'z' - :param numpy.ndarray moment: The vector dipole moment (vertical) - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location - """ - - if isinstance(orientation, str): - assert orientation.upper() in ['X', 'Y', 'Z'], ( - "orientation must be 'x', 'y', or 'z' or a vector not {}" - .format(orientation) - ) - elif (not np.allclose(np.r_[1., 0., 0.], orientation) or - not np.allclose(np.r_[0., 1., 0.], orientation) or - not np.allclose(np.r_[0., 0., 1.], orientation)): - warnings.warn( - 'Arbitrary trasnmitter orientations ({}) not thouroughly tested ' - 'Pull request on a test anyone? bueller?'.format(orientation) - ) - - if isinstance(component, str): - assert component.upper() in ['X', 'Y', 'Z'], ( - "component must be 'x', 'y', or 'z' or a vector not {}" - .format(component) - ) - elif (not np.allclose(np.r_[1., 0., 0.], component) or - not np.allclose(np.r_[0., 1., 0.], component) or - not np.allclose(np.r_[0., 0., 1.], component)): - warnings.warn( - 'Arbitrary receiver orientations ({}) not thouroughly tested ' - 'Pull request on a test anyone? bueller?' - ).format(component) - - if isinstance(orientation, str): - orientation = orientationDict[orientation.upper()] - - if isinstance(component, str): - component = orientationDict[component.upper()] - - assert np.linalg.norm(orientation, 2) == 1., ( - "orientation must be a unit vector. " - "Use 'moment=X to scale source fields" - ) - - if np.linalg.norm(component, 2) != 1.: - warnings.warn( - 'The magnitude of the receiver component vector is > 1, ' - ' it is {}. The receiver fields will be scaled.' - .format(np.linalg.norm(component, 2)) - ) - - srcLoc = np.atleast_2d(srcLoc) - component = np.atleast_2d(component) - obsLoc = np.atleast_2d(obsLoc) - orientation = np.atleast_2d(orientation) - - nObs = obsLoc.shape[0] - nSrc = int(srcLoc.size / 3.) - - # use outer product to construct an array of [x_src, y_src, z_src] - - m = moment*orientation.repeat(nObs, axis=0) - B = [] - - for i in range(nSrc): - srcLoc = srcLoc[i, np.newaxis].repeat(nObs, axis=0) - rx = component.repeat(nObs, axis=0) - dR = obsLoc - srcLoc - r = np.sqrt((dR**2).sum(axis=1)) - - # mult each element and sum along the axis (vector dot product) - m_dot_dR_div_r2 = (m * dR).sum(axis=1) / (r**2) - - # multiply the scalar m_dot_dR by the 3D vector r - rvec_m_dot_dR_div_r2 = np.vstack([m_dot_dR_div_r2 * dR[:, i] for - i in range(3)]).T - inside = (3. * rvec_m_dot_dR_div_r2) - m - - # dot product with rx orientation - inside_dot_rx = (inside * rx).sum(axis=1) - front = (mu/(4.* pi * r**3)) - - B.append(Utils.mkvc(front * inside_dot_rx)) - - return np.vstack(B).T diff --git a/geoana/spatial.py b/geoana/spatial.py index d6d5f99c..07d429db 100644 --- a/geoana/spatial.py +++ b/geoana/spatial.py @@ -43,6 +43,18 @@ def vector_dot(xyz, vector): :param numpy.array vector: vector (1 x 3) """ assert len(vector) == 3, "vector should be length 3" - return np.hstack([ - vector[0]*xyz[:, 0], vector[1]*xyz[:, 1], vector[2]*xyz[:, 2] - ]) + return vector[0]*xyz[:, 0] + vector[1]*xyz[:, 1] + vector[2]*xyz[:, 2] + + +def repeat_scalar(scalar, dim=3): + """ + Repeat a spatially distributed scalar value dim times to simplify + multiplication with a vector. + """ + assert len(scalar) in scalar.shape, ( + "input must be a scalar. The shape you provided is {}".format( + scalar.shape + ) + ) + + return np.kron(np.ones((1, dim)), np.atleast_2d(scalar).T) diff --git a/requirements_dev.txt b/requirements_dev.txt new file mode 100644 index 00000000..5cc9a934 --- /dev/null +++ b/requirements_dev.txt @@ -0,0 +1,2 @@ +discretize +sphinx_rtd_theme diff --git a/tests/test_em_fdem.py b/tests/test_em_fdem.py deleted file mode 100644 index 987fc862..00000000 --- a/tests/test_em_fdem.py +++ /dev/null @@ -1,73 +0,0 @@ -from __future__ import absolute_import -from __future__ import division -from __future__ import print_function -from __future__ import unicode_literals - -import unittest - -from geoana.em import fdem - - -def E_from_EDWS(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1., t=0.): - """E_from_EDWS - Computing the analytic electric fields (E) from an electrical dipole in a wholespace - - You have the option of computing E for multiple frequencies at a single reciever location - or a single frequency at multiple locations - - :param numpy.array XYZ: reciever locations at which to evaluate E - :param float epsr: relative permitivitty value (unitless), default is 1.0 - :rtype: numpy.array - :return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and frequencies. - """ - - mu = mu_0*(1+kappa) - epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omega(f)*epsilon - - XYZ = Utils.asArray_N_x_Dim(XYZ, 3) - # Check - if XYZ.shape[0] > 1 & f.shape[0] > 1: - raise Exception('I/O type error: For multiple field locations only a single frequency can be specified.') - - dx = XYZ[:, 0] - srcLoc[0] - dy = XYZ[:, 1] - srcLoc[1] - dz = XYZ[:, 2] - srcLoc[2] - - r = np.sqrt(dx**2. + dy**2. + dz**2.) - # k = np.sqrt( -1j*2.*pi*f*mu*sig ) - k = np.sqrt(omega(f)**2. * mu * epsilon - 1j * omega(f) * mu * sig) - - front = current * length / (4.*pi*sig_hat* r**3) * np.exp(-1j*k*r) - mid = -k**2 * r**2 + 3*1j*k*r + 3 - - if orientation.upper() == 'X': - Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ey = front*(dx*dy / r**2)*mid - Ez = front*(dx*dz / r**2)*mid - return Ex, Ey, Ez - - elif orientation.upper() == 'Y': - # x--> y, y--> z, z-->x - Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ez = front*(dy*dz / r**2)*mid - Ex = front*(dy*dx / r**2)*mid - return Ex, Ey, Ez - - elif orientation.upper() == 'Z': - # x --> z, y --> x, z --> y - Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) - Ex = front*(dz*dx / r**2)*mid - Ey = front*(dz*dy / r**2)*mid - return Ex, Ey, Ez - - -class TestFDEM(unittest.TestCase): - - def test_vector(self): - edws = fdem.ElectricDipole_WholeSpace() - - - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_em_fdem_electric_dipole.py b/tests/test_em_fdem_electric_dipole.py new file mode 100644 index 00000000..d016574b --- /dev/null +++ b/tests/test_em_fdem_electric_dipole.py @@ -0,0 +1,187 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + +import unittest +import numpy as np + +from scipy.constants import mu_0, epsilon_0 +from geoana.em import fdem +from discretize.utils import ndgrid, asArray_N_x_Dim + + +def E_from_EDWS( + XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., + epsr=1., t=0. +): + """E_from_EDWS + Computing the analytic electric fields (E) from an electrical dipole in + a wholespace + - You have the option of computing E for multiple frequencies at a single + reciever location + or a single frequency at multiple locations + + :param numpy.array XYZ: reciever locations at which to evaluate E + :param float epsr: relative permitivitty value (unitless), default is 1.0 + :rtype: numpy.array + :return: Ex, Ey, Ez: arrays containing all 3 components of E evaluated at the specified locations and frequencies. + """ + + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*fdem.omega(f)*epsilon + + XYZ = asArray_N_x_Dim(XYZ, 3) + + dx = XYZ[:, 0] - srcLoc[0] + dy = XYZ[:, 1] - srcLoc[1] + dz = XYZ[:, 2] - srcLoc[2] + + r = np.sqrt(dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*pi*f*mu*sig ) + k = np.sqrt( + fdem.omega(f)**2. * mu * epsilon - 1j * fdem.omega(f) * mu * sig + ) + + front = current * length / (4.*np.pi*sig_hat * r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r-1.)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r-1.)) + Ez = front*(dy*dz / r**2)*mid + Ex = front*(dy*dx / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + + +class TestFDEMdipole(unittest.TestCase): + + def test_defaults(self): + edws = fdem.ElectricDipoleWholeSpace() + assert(edws.sigma == 1) + assert(edws.mu == mu_0) + assert(edws.epsilon == epsilon_0) + assert(np.all(edws.orientation == np.r_[1., 0., 0.])) + assert(edws.length == 1.) + assert(np.all(edws.location == np.r_[0., 0., 0.])) + assert(edws.frequency == 1.) + + def compare_fields(name, field, ftest): + + def check_component(name, f, ftest): + geoana_norm = np.linalg.norm(f) + test_norm = np.linalg.norm(ftest) + diff = np.linalg.norm(f-ftest) + passed = np.allclose(f, ftest) + print( + "Testing {} ... geoana: {:1.4e}, compare: {:1.4e}, " + "diff: {:1.4e}, passed?: {}".format( + name, geoana_norm, test_norm, diff, passed + ) + ) + return passed + + for i, orientation in enumerate(['x', 'y', 'z']): + for component in ['real', 'imag']: + check_component( + orientation + '_' + component, + getattr(field[:, i], component), + getattr(ftest[:, i], component) + ) + + def electric_dipole_e(self, orientation): + sigma = np.exp(np.random.randn(1)) + frequency = np.random.rand(1)*1e6 + edws = fdem.ElectricDipoleWholeSpace( + orientation=orientation, + sigma=sigma, + frequency=frequency + ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + xyz = ndgrid([x, y, z]) + + extest, eytest, eztest = E_from_EDWS( + xyz, edws.location, edws.sigma, edws.frequency, + orientation=orientation.upper() + ) + + e = edws.electric_field(xyz) + print( + "\n\nTesting Electric Dipole {} orientation\n".format(orientation) + ) + + self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + + def test_electric_dipole_x_e(self): + self.electric_dipole_e("x") + + def test_electric_dipole_y_e(self): + self.electric_dipole_e("y") + + def test_electric_dipole_z_e(self): + self.electric_dipole_e("z") + + def test_electric_dipole_tilted_e(self): + + orientation = np.random.rand(3) + orientation = orientation / np.linalg.norm(orientation) + + edws = fdem.ElectricDipoleWholeSpace( + orientation=orientation + ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + + xyz = ndgrid([x, y, z]) + + extest0, eytest0, eztest0 = E_from_EDWS( + xyz, edws.location, edws.sigma, edws.frequency, + orientation='X' + ) + extest1, eytest1, eztest1 = E_from_EDWS( + xyz, edws.location, edws.sigma, edws.frequency, + orientation='Y' + ) + extest2, eytest2, eztest2 = E_from_EDWS( + xyz, edws.location, edws.sigma, edws.frequency, + orientation='Z' + ) + + extest = ( + orientation[0]*extest0 + orientation[1]*extest1 + orientation[2]*extest2 + ) + eytest = ( + orientation[0]*eytest0 + orientation[1]*eytest1 + orientation[2]*eytest2 + ) + eztest = ( + orientation[0]*eztest0 + orientation[1]*eztest1 + orientation[2]*eztest2 + ) + + e = edws.electric_field(xyz) + print( + "\n\nTesting Electric Dipole {} orientation\n".format("45 degree") + ) + + self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_em_fdem_magnetic_dipole.py b/tests/test_em_fdem_magnetic_dipole.py new file mode 100644 index 00000000..c7a48181 --- /dev/null +++ b/tests/test_em_fdem_magnetic_dipole.py @@ -0,0 +1,112 @@ +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals + +import unittest +import numpy as np + +from scipy.constants import mu_0, epsilon_0 +from geoana.em import fdem +from discretize.utils import ndgrid, asArray_N_x_Dim + + +def MagneticDipoleFields(srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0): + """ + Calculate the vector potential of a set of magnetic dipoles + at given locations 'ref. ' + + .. math:: + + B = \frac{\mu_0}{4 \pi r^3} \left( \frac{3 \vec{r} (\vec{m} \cdot + \vec{r})}{r^2}) + - \vec{m} + \right) \cdot{\hat{rx}} + + :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) + :param numpy.ndarray obsLoc: Where the potentials will be calculated + (x, y, z) + :param str component: The component to calculate - 'x', 'y', or 'z' + :param numpy.ndarray moment: The vector dipole moment (vertical) + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if isinstance(orientation, str): + assert orientation.upper() in ['X', 'Y', 'Z'], ( + "orientation must be 'x', 'y', or 'z' or a vector not {}" + .format(orientation) + ) + elif (not np.allclose(np.r_[1., 0., 0.], orientation) or + not np.allclose(np.r_[0., 1., 0.], orientation) or + not np.allclose(np.r_[0., 0., 1.], orientation)): + warnings.warn( + 'Arbitrary trasnmitter orientations ({}) not thouroughly tested ' + 'Pull request on a test anyone? bueller?'.format(orientation) + ) + + if isinstance(component, str): + assert component.upper() in ['X', 'Y', 'Z'], ( + "component must be 'x', 'y', or 'z' or a vector not {}" + .format(component) + ) + elif (not np.allclose(np.r_[1., 0., 0.], component) or + not np.allclose(np.r_[0., 1., 0.], component) or + not np.allclose(np.r_[0., 0., 1.], component)): + warnings.warn( + 'Arbitrary receiver orientations ({}) not thouroughly tested ' + 'Pull request on a test anyone? bueller?' + ).format(component) + + if isinstance(orientation, str): + orientation = orientationDict[orientation.upper()] + + if isinstance(component, str): + component = orientationDict[component.upper()] + + assert np.linalg.norm(orientation, 2) == 1., ( + "orientation must be a unit vector. " + "Use 'moment=X to scale source fields" + ) + + if np.linalg.norm(component, 2) != 1.: + warnings.warn( + 'The magnitude of the receiver component vector is > 1, ' + ' it is {}. The receiver fields will be scaled.' + .format(np.linalg.norm(component, 2)) + ) + + srcLoc = np.atleast_2d(srcLoc) + component = np.atleast_2d(component) + obsLoc = np.atleast_2d(obsLoc) + orientation = np.atleast_2d(orientation) + + nObs = obsLoc.shape[0] + nSrc = int(srcLoc.size / 3.) + + # use outer product to construct an array of [x_src, y_src, z_src] + + m = moment*orientation.repeat(nObs, axis=0) + B = [] + + for i in range(nSrc): + srcLoc = srcLoc[i, np.newaxis].repeat(nObs, axis=0) + rx = component.repeat(nObs, axis=0) + dR = obsLoc - srcLoc + r = np.sqrt((dR**2).sum(axis=1)) + + # mult each element and sum along the axis (vector dot product) + m_dot_dR_div_r2 = (m * dR).sum(axis=1) / (r**2) + + # multiply the scalar m_dot_dR by the 3D vector r + rvec_m_dot_dR_div_r2 = np.vstack([m_dot_dR_div_r2 * dR[:, i] for + i in range(3)]).T + inside = (3. * rvec_m_dot_dR_div_r2) - m + + # dot product with rx orientation + inside_dot_rx = (inside * rx).sum(axis=1) + front = (mu/(4.* pi * r**3)) + + B.append(Utils.mkvc(front * inside_dot_rx)) + + return np.vstack(B).T From 6fa7798fce28907d59816f57cc4d739ee9358c34 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 26 Jun 2017 15:15:20 -0700 Subject: [PATCH 09/41] Add SimPEG to dev requirements --- requirements_dev.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements_dev.txt b/requirements_dev.txt index 5cc9a934..b7b10cd8 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,2 +1,3 @@ discretize sphinx_rtd_theme +SimPEG From de22e01470ba81d9507964e51299e05d0ca483ac Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 26 Jun 2017 20:30:37 -0700 Subject: [PATCH 10/41] test against SimPEG --- geoana/em/base.py | 2 +- geoana/em/fdem.py | 4 +- tests/test_em_fdem_electric_dipole.py | 165 ++++++++++++++++++++++++-- 3 files changed, 159 insertions(+), 12 deletions(-) diff --git a/geoana/em/base.py b/geoana/em/base.py index 231a2872..b6655acd 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -108,7 +108,7 @@ def cross_orientation(self, xyz): Take the cross product between a grid and the orientation of the dipole """ orientation = np.kron( - np.atleast_2d(self.orientation), np.ones(xyz.shape[0]).T + np.atleast_2d(self.orientation), np.ones((xyz.shape[0], 1)) ) return np.cross(xyz, orientation) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 7df35295..6f57b44d 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -238,14 +238,14 @@ def magnetic_field(self, xyz): """ dxyz = self.vector_distance(xyz) r = self.distance(xyz) - r = np.kron(np.ones(1, 3), np.atleast_2d(r).T) + r = spatial.repeat_scalar(r) kr = self.wave_number * r front = ( self.current * self.length / (4 * np.pi * r**2) * (1j * kr + 1) * np.exp(-1j * kr) ) - return front * self.cross_orientation(xyz) / r + return - front * self.cross_orientation(xyz) / r def magnetic_flux_density(self, xyz): """ diff --git a/tests/test_em_fdem_electric_dipole.py b/tests/test_em_fdem_electric_dipole.py index d016574b..f914fecf 100644 --- a/tests/test_em_fdem_electric_dipole.py +++ b/tests/test_em_fdem_electric_dipole.py @@ -8,7 +8,10 @@ from scipy.constants import mu_0, epsilon_0 from geoana.em import fdem -from discretize.utils import ndgrid, asArray_N_x_Dim +import discretize + +from SimPEG.EM import FDEM +from SimPEG import Maps def E_from_EDWS( @@ -32,7 +35,7 @@ def E_from_EDWS( epsilon = epsilon_0*epsr sig_hat = sig + 1j*fdem.omega(f)*epsilon - XYZ = asArray_N_x_Dim(XYZ, 3) + XYZ = discretize.utils.asArray_N_x_Dim(XYZ, 3) dx = XYZ[:, 0] - srcLoc[0] dy = XYZ[:, 1] - srcLoc[1] @@ -95,17 +98,19 @@ def check_component(name, f, ftest): ) return passed + passed = [] for i, orientation in enumerate(['x', 'y', 'z']): for component in ['real', 'imag']: - check_component( + passed.append(check_component( orientation + '_' + component, getattr(field[:, i], component), getattr(ftest[:, i], component) - ) + )) + return all(passed) def electric_dipole_e(self, orientation): - sigma = np.exp(np.random.randn(1)) - frequency = np.random.rand(1)*1e6 + sigma = np.random.random_integers(1) + frequency = np.random.random_integers(1) edws = fdem.ElectricDipoleWholeSpace( orientation=orientation, sigma=sigma, @@ -114,7 +119,7 @@ def electric_dipole_e(self, orientation): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = ndgrid([x, y, z]) + xyz = discretize.utils.ndgrid([x, y, z]) extest, eytest, eztest = E_from_EDWS( xyz, edws.location, edws.sigma, edws.frequency, @@ -126,7 +131,8 @@ def electric_dipole_e(self, orientation): "\n\nTesting Electric Dipole {} orientation\n".format(orientation) ) - self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + passed = self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + self.assertTrue(passed) def test_electric_dipole_x_e(self): self.electric_dipole_e("x") @@ -149,7 +155,7 @@ def test_electric_dipole_tilted_e(self): y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = ndgrid([x, y, z]) + xyz = discretize.utils.ndgrid([x, y, z]) extest0, eytest0, eztest0 = E_from_EDWS( xyz, edws.location, edws.sigma, edws.frequency, @@ -182,6 +188,147 @@ def test_electric_dipole_tilted_e(self): self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) +class TestFDEMdipole_SimPEG(unittest.TestCase): + + tol = 1e-1 # error must be an order of magnitude smaller than results + + # put the source at the center + + def getFaceSrc(self, mesh): + csx = mesh.hx.min() + csz = mesh.hz.min() + srcInd = ( + (mesh.gridFz[:, 0] < csx) & + (mesh.gridFz[:, 2] < csz/2.) & (mesh.gridFz[:, 2] > -csz/2.) + ) + + src_vecz = np.zeros(mesh.nFz, dtype=complex) + src_vecz[srcInd] = 1. + + return np.hstack( + [np.zeros(mesh.vnF[:2].sum(), dtype=complex), src_vecz] + ) + + def getProjections(self, mesh): + ignore_inside_radius = 5*mesh.hx.min() + ignore_outside_radius = 40*mesh.hx.min() + + def ignoredGridLocs(grid): + return ( + ( + grid[:, 0]**2 + grid[:, 1]**2 + grid[:, 2]**2 < + ignore_inside_radius**2 + ) | ( + grid[:, 0]**2 + grid[:, 1]**2 + grid[:, 2]**2 > + ignore_outside_radius**2 + ) + ) + + # Faces + ignore_me_Fx = ignoredGridLocs(mesh.gridFx) + ignore_me_Fz = ignoredGridLocs(mesh.gridFz) + ignore_me = np.hstack([ignore_me_Fx, ignore_me_Fz]) + keep_me = np.array(~ignore_me, dtype=float) + Pf = discretize.utils.sdiag(keep_me) + + # Edges + ignore_me_Ey = ignoredGridLocs(mesh.gridEy) + keep_me_Ey = np.array(~ignore_me_Ey, dtype=float) + Pe = discretize.utils.sdiag(keep_me_Ey) + + return Pf, Pe + + def test_e_dipole_v_SimPEG(self): + + def compare_w_SimPEG(name, geoana, simpeg): + + norm_geoana = np.linalg.norm(geoana) + norm_simpeg = np.linalg.norm(simpeg) + diff = np.linalg.norm(geoana - simpeg) + passed = diff < self.tol * 0.5 * (norm_geoana + norm_simpeg) + print( + " {} ... geoana: {:1.4e}, SimPEG: {:1.4e}, diff: {:1.4e}, " + "passed?: {}".format( + name, norm_geoana, norm_simpeg, diff, passed + ) + ) + + return passed + + print("\n\nComparing Electric dipole with SimPEG") + + sigma_back = 1e-1 + freqs = np.r_[10., 100.] + + csx, ncx, npadx = 0.5, 40, 20 + ncy = 1 + csz, ncz, npadz = 0.5, 40, 20 + + hx = discretize.utils.meshTensor( + [(csx, ncx), (csx, npadx, 1.2)] + ) + hy = 2*np.pi / ncy * np.ones(ncy) + hz = discretize.utils.meshTensor( + [(csz, npadz, -1.2), (csz, ncz), (csz, npadz, 1.2)] + ) + + mesh = discretize.CylMesh([hx, hy, hz], x0='00C') + + s_e = self.getFaceSrc(mesh) + prob = FDEM.Problem3D_h(mesh, sigmaMap=Maps.IdentityMap(mesh)) + srcList = [FDEM.Src.RawVec_e([], f, s_e) for f in freqs] + survey = FDEM.Survey(srcList) + + prob.pair(survey) + + fields = prob.fields(sigma_back*np.ones(mesh.nC)) + + length = mesh.hz.min() + current = np.pi * csx**2 + + edws = fdem.ElectricDipoleWholeSpace( + sigma=sigma_back, length=length, current=current, orientation="z" + ) + + Pf, Pe = self.getProjections(mesh) + + j_passed = [] + h_passed = [] + for i, f in enumerate(freqs): + edws.frequency = f + + j_xz = [] + for j, component in zip([0, 2], ['x', 'z']): + grid = getattr(mesh, "gridF{}".format(component)) + j_xz.append( + edws.current_density(grid)[:, j + ] + ) + j_geoana = np.hstack(j_xz) + h_geoana = edws.magnetic_field(mesh.gridEy)[:, 1] + + P_j_geoana = Pf*j_geoana + P_j_simpeg = Pf*discretize.utils.mkvc(fields[srcList[i], 'j']) + + P_h_geoana = Pe*h_geoana + P_h_simpeg = Pe*discretize.utils.mkvc(fields[srcList[i], 'h']) + + print("Testing {} Hz".format(f)) + + for comp in ['real', 'imag']: + j_passed.append(compare_w_SimPEG( + 'J {}'.format(comp), + getattr(P_j_geoana, comp), + getattr(P_j_simpeg, comp) + )) + h_passed.append(compare_w_SimPEG( + 'H {}'.format(comp), + getattr(P_h_geoana, comp), + getattr(P_h_simpeg, comp) + )) + assert(all(j_passed)) + assert(all(h_passed)) + if __name__ == '__main__': unittest.main() From ab94986e427a1574f82c4567ee473b1a1506b676 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 16:33:47 -0700 Subject: [PATCH 11/41] properties upgrades in the earthquake code --- geoana/earthquake/oksar.py | 16 ++++++++++++++-- tests/test_earthquake_oksar.py | 3 +-- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/geoana/earthquake/oksar.py b/geoana/earthquake/oksar.py index 11451e7e..89ed6476 100644 --- a/geoana/earthquake/oksar.py +++ b/geoana/earthquake/oksar.py @@ -23,7 +23,17 @@ import matplotlib.pyplot as plt -class EarthquakeInterferogram(properties.UidModel): +class EarthquakeInterferogram(properties.HasProperties): + title = properties.String( + 'name of the earthquake', + required=True + ) + + description = properties.String( + 'description of the event', + required=False + ) + location = properties.Vector2( 'interferogram location (bottom N, left E)', required=True @@ -229,6 +239,8 @@ def get_LOS_vector(self, locations): calculate beta - the angle at earth center between reference point and satellite nadir """ + if not isinstance(locations, list): + locations = [locations] utmZone = self.location_UTM_zone refPoint = vmath.Vector3(self.ref.x, self.ref.y, 0) @@ -283,7 +295,7 @@ def get_LOS_vector(self, locations): los_y = -np.sin(satAzimuth * DEG2RAD) * np.cos(satIncidence * DEG2RAD) los_z = np.sin(satIncidence * DEG2RAD) - return vmath.Vector3(los_x, los_y, los_z) + return vmath.Vector3Array([los_x, los_y, los_z]) @staticmethod def _ang_to_gc(x, y, origx, origy, satAzimuth): diff --git a/tests/test_earthquake_oksar.py b/tests/test_earthquake_oksar.py index be451f3e..a5c97a12 100644 --- a/tests/test_earthquake_oksar.py +++ b/tests/test_earthquake_oksar.py @@ -21,9 +21,8 @@ def test_los(self): dinar.satellite_azimuth = 192 dinar.location_UTM_zone = 35 - utmLoc = vmath.Vector3([706216.0606], [4269238.9999], [0]) + utmLoc = vmath.Vector3([706216.0606, 4269238.9999, 0]) # refPoint = vmath.Vector3(dinar.ref.x, dinar.ref.y, 0) - LOS = dinar.get_LOS_vector(utmLoc) # compare against fortran code. From 94415f9049c7834ba36ebde17f2c9bbb125e3cf5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 19:06:58 -0700 Subject: [PATCH 12/41] add earthquakes to the init, add nightly builds of python to travis testing --- .travis.yml | 6 +++++- geoana/__init__.py | 2 ++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index a381b053..ee9d0f4b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,7 +1,11 @@ language: python python: - - 3.4 + - 3.6 - 2.7 + - "nightly" + +allowed_failures: + - "nightly" sudo: false diff --git a/geoana/__init__.py b/geoana/__init__.py index 1edc33ac..fc23e4d4 100644 --- a/geoana/__init__.py +++ b/geoana/__init__.py @@ -1 +1,3 @@ +from . import earthquake from . import em + From 85cd421ec4684d9a615e556d525677c82bfcf79e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 19:17:41 -0700 Subject: [PATCH 13/41] fix up the readme anchor --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 76ad6ca4..4698ac22 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,4 @@ -.. _index: +.. _geoana: geoana ====== From fcd24cd405380d66c458387a52b459ff8c4bd5ef Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 19:26:43 -0700 Subject: [PATCH 14/41] Add testing on appveyor --- appveyor.yml | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 appveyor.yml diff --git a/appveyor.yml b/appveyor.yml new file mode 100644 index 00000000..e69de29b From eb22df72115e2dcd464146c9a66317f3e29b3e06 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 19:30:02 -0700 Subject: [PATCH 15/41] cleanup appveyor.yml --- appveyor.yml | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/appveyor.yml b/appveyor.yml index e69de29b..9906bc68 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -0,0 +1,45 @@ +# AppVeyor.com is a Continuous Integration service to build and run tests under Windows + +build: off + +environment: + matrix: + - PYTHON: 2.7 + MINICONDA: C:\Miniconda + PYTHON_ARCH: 32 + + - PYTHON: 3.5 + MINICONDA: C:\Miniconda3 + PYTHON_ARCH: 32 + + - PYTHON: 3.6 + MINICONDA: C:\Miniconda3 + PYTHON_ARCH: 32 + + - PYTHON: 2.7 + MINICONDA: C:\Miniconda + PYTHON_ARCH: 64 + + - PYTHON: 3.5 + MINICONDA: C:\Miniconda3 + PYTHON_ARCH: 64 + + - PYTHON: 3.6 + MINICONDA: C:\Miniconda3 + PYTHON_ARCH: 64 + +init: + - "ECHO %PYTHON% %MINICONDA%" + +install: + - "set PATH=%MINICONDA%;%MINICONDA%\\Scripts;%PATH%" + - conda config --set always_yes yes --set changeps1 no + - conda update -q conda + - conda info -a + - "conda create -q -n test-environment python=%PYTHON% numpy scipy matplotlib cython ipython pillow wheel" + - activate test-environment + - pip install -r requirements_dev.txt + - python setup.py install + +test_script: + - nosetests tests -v -s From f29635f95439f11ac0a0ebf71d4cbe3dbc7169c8 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 19:40:14 -0700 Subject: [PATCH 16/41] boilerplate stuff --- .bumpversion.cfg | 4 ++++ requirements.txt | 9 +------- setup.py | 57 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 62 insertions(+), 8 deletions(-) create mode 100644 .bumpversion.cfg create mode 100644 setup.py diff --git a/.bumpversion.cfg b/.bumpversion.cfg new file mode 100644 index 00000000..bb809cd0 --- /dev/null +++ b/.bumpversion.cfg @@ -0,0 +1,4 @@ +[bumpversion] +current_version = 0.0.1 +files = setup.py geoana/__init__.py docs/conf.py + diff --git a/requirements.txt b/requirements.txt index 56430129..d6e1198b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1 @@ -numpy -scipy -matplotlib -sphinx -ipywidgets -requests -utm -properties +-e . diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..5c30b42c --- /dev/null +++ b/setup.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python +from __future__ import print_function +"""geoana + +Interactive geoscience (mostly) analytic functions. +""" + +from distutils.core import setup +from setuptools import find_packages + + +CLASSIFIERS = [ + 'Development Status :: 4 - Beta', + 'Intended Audience :: Education', + 'Intended Audience :: Science/Research', + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python', + 'Topic :: Scientific/Engineering', + 'Topic :: Scientific/Engineering :: Mathematics', + 'Topic :: Scientific/Engineering :: Physics', + 'Operating System :: Microsoft :: Windows', + 'Operating System :: POSIX', + 'Operating System :: Unix', + 'Operating System :: MacOS', + 'Natural Language :: English', +] + +with open('README.rst') as f: + LONG_DESCRIPTION = ''.join(f.readlines()) + +setup( + name = 'geoana', + version = '0.0.1', + packages = find_packages(), + install_requires = [ + 'future', + 'numpy>=1.7', + 'scipy>=0.13', + 'matplotlib', + 'requests', + 'ipywidgets', + 'properties[math]', + 'jupyter', + 'utm' + ], + author = 'SimPEG developers', + author_email = 'lindseyheagy@gmail.com', + description = 'geoana', + long_description = LONG_DESCRIPTION, + keywords = 'geophysics, electromagnetics', + url = 'http://simpeg.xyz', + download_url = 'https://github.com/simpeg/geoana', + classifiers=CLASSIFIERS, + platforms = ['Windows', 'Linux', 'Solaris', 'Mac OS-X', 'Unix'], + license='MIT License', + use_2to3 = False, +) From 4f5558f6777ae1969b14247011bef09cda6fe002 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 21:30:59 -0700 Subject: [PATCH 17/41] suppress non local image warnings, makefile for tests, clean, etc --- Makefile | 34 ++++++++++++++++++++++++++++++++++ docs/conf.py | 14 ++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 Makefile diff --git a/Makefile b/Makefile new file mode 100644 index 00000000..ba895dc5 --- /dev/null +++ b/Makefile @@ -0,0 +1,34 @@ +PACKAGE_NAME=geoana + +.PHONY: install coverage lint graphs tests docs clean deploy + +install: + python setup.py install + +coverage: + nosetests --logging-level=INFO --with-coverage --cover-package=discretize --cover-html + open cover/index.html + +lint: + pylint $(PACKAGE_NAME) + +lint-html: + pylint --output-format=html $(PACKAGE_NAME) > pylint.html + +graphs: + pyreverse -my -A -o pdf -p discretize discretize/**.py discretize/**/**.py + +tests: + nosetests --logging-level=INFO + +docs: + cd docs;make html + +clean_pyc: + find . -name "*.pyc" | xargs -I {} rm -v "{}" + +clean: clean_pyc + cd docs;make clean + +deploy: + python setup.py sdist upload diff --git a/docs/conf.py b/docs/conf.py index cbe511af..6acf383f 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -299,3 +299,17 @@ # Example configuration for intersphinx: refer to the Python standard library. intersphinx_mapping = {'https://docs.python.org/': None} + +def supress_nonlocal_image_warn(): + import sphinx.environment + sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn + + +def _supress_nonlocal_image_warn(self, msg, node, **kwargs): + from docutils.utils import get_source_line + + if not msg.startswith('nonlocal image URI found:'): + self._warnfunc(msg, '{0!s}:{1!s}'.format(*get_source_line(node))) + +supress_nonlocal_image_warn() + From d0654ecc02fa56352377f2d1044b738f7b43092f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 27 Jun 2017 21:42:57 -0700 Subject: [PATCH 18/41] supress non local image warning in sphinx, cleanup docstrings in fdem --- docs/conf.py | 1 + geoana/__init__.py | 5 +++++ geoana/em/base.py | 2 +- geoana/em/fdem.py | 36 +++++++++++++++++------------------- 4 files changed, 24 insertions(+), 20 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 6acf383f..765dab16 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -300,6 +300,7 @@ # Example configuration for intersphinx: refer to the Python standard library. intersphinx_mapping = {'https://docs.python.org/': None} + def supress_nonlocal_image_warn(): import sphinx.environment sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn diff --git a/geoana/__init__.py b/geoana/__init__.py index fc23e4d4..ae43d999 100644 --- a/geoana/__init__.py +++ b/geoana/__init__.py @@ -1,3 +1,8 @@ from . import earthquake from . import em +__version__ = '0.0.1' +__author__ = 'SimPEG developers' +__license__ = 'MIT' +__copyright__ = 'Copyright 2017 SimPEG developers' + diff --git a/geoana/em/base.py b/geoana/em/base.py index b6655acd..2871564e 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -141,7 +141,7 @@ class BaseElectricDipole(BaseDipole): ) current = properties.Float( - "size of the injected current (A)", + "magnitude of the injected current (A)", default=1.0, min=0.0 ) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 6f57b44d..768eafe4 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -138,19 +138,20 @@ def sigma_hat(self): @property def wave_number(self): """ - Wavenumber of an electromagnetic wave in a medium with constant physical - properties + Wavenumber of an electromagnetic wave in a medium with constant + physical properties .. math:: - k = \sqrt{\omega**2 \mu \varepsilon - i \omega \mu \sigma} """ + k = \sqrt{\omega**2 \mu \varepsilon - i \omega \mu \sigma} + """ return wave_number(self.frequency, self.sigma, self.mu) @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 + Distance at which an em wave has decayed by a factor of :math:`1/e` in + a medium with constant physical properties .. math:: @@ -159,16 +160,7 @@ def skin_depth(self): return skin_depth(self.frequency, self.sigma, self.mu) -class BaseFDEMDipoleWholeSpace(BaseFDEM): - """ - Base FDEM Dipole - """ - pass - - -class ElectricDipoleWholeSpace( - BaseElectricDipole, BaseFDEMDipoleWholeSpace -): +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 @@ -177,7 +169,9 @@ class ElectricDipoleWholeSpace( .. math:: - \mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r} - \mathbf{r}_s)\mathbf{\hat{u}} + \mathbf{J}(\mathbf{r}) = I ds \delta(\mathbf{r} + - \mathbf{r}_s)\mathbf{\hat{u}} + """ def vector_potential(self, xyz): """ @@ -202,7 +196,8 @@ def electric_field(self, xyz): .. math:: - \mathbf{E} = \frac{1}{\hat{\sigma}} \nabla \nabla \cdot \mathbf{A} - i \omega \mu \mathbf{A} + \mathbf{E} = \frac{1}{\hat{\sigma}} \nabla \nabla \cdot \mathbf{A} + - i \omega \mu \mathbf{A} """ dxyz = self.vector_distance(xyz) @@ -211,7 +206,7 @@ def electric_field(self, xyz): kr = self.wave_number * r front = ( - (self.current * self.length) / (4 * np.pi * self.sigma * r**3 ) * + (self.current * self.length) / (4 * np.pi * self.sigma * r**3) * np.exp(-1j * kr) ) symmetric_term = ( @@ -222,9 +217,12 @@ def electric_field(self, xyz): (kr**2 - 1j*kr - 1) * np.kron(self.orientation, np.ones((dxyz.shape[0], 1))) ) - return front * ( symmetric_term + oriented_term ) + return front * (symmetric_term + oriented_term) def current_density(self, xyz): + """ + Current density due to a harmonic electric dipole + """ return self.sigma * self.electric_field(xyz) def magnetic_field(self, xyz): From 330dd8747bffa2bd6ad6890eb08d3abcbedb39e3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 28 Jun 2017 12:45:10 -0700 Subject: [PATCH 19/41] add sphinx to dev requirements --- requirements_dev.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements_dev.txt b/requirements_dev.txt index b7b10cd8..23d1fb05 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,3 +1,4 @@ discretize sphinx_rtd_theme SimPEG +sphinx From 9464f98766b0f01eede90550e1893ff58739ba42 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 28 Jun 2017 12:52:43 -0700 Subject: [PATCH 20/41] fix nightly build to be an allowed failure --- .travis.yml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index ee9d0f4b..b6489bcb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,8 +4,9 @@ python: - 2.7 - "nightly" -allowed_failures: - - "nightly" +matrix: + allow_failures: + - "nightly" sudo: false From 2668a3dcd7da124d9b52dd76fbe5e0dbbf657f44 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 30 Jun 2017 08:41:01 -0700 Subject: [PATCH 21/41] start of TDEM electric dipole --- .travis.yml | 5 -- geoana/em/base.py | 42 ---------- geoana/em/fdem.py | 15 +++- geoana/em/tdem.py | 206 ++++++++++++++++++++++++++++++++++------------ 4 files changed, 166 insertions(+), 102 deletions(-) diff --git a/.travis.yml b/.travis.yml index b6489bcb..9aff16c8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,11 +2,6 @@ language: python python: - 3.6 - 2.7 - - "nightly" - -matrix: - allow_failures: - - "nightly" sudo: false diff --git a/geoana/em/base.py b/geoana/em/base.py index 2871564e..ca7e2c50 100644 --- a/geoana/em/base.py +++ b/geoana/em/base.py @@ -10,32 +10,6 @@ from .. import spatial -def peak_time(z, sigma, mu=mu_0): - """ - Time at which the maximum signal amplitude is observed at a particular - location for a transient plane wave through a homogeneous medium. - - See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html - - :param z float: distance from source (m) - :param sigma float: electrical conductivity (S/m) - :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m - """ - return (mu * sigma * z**2)/6. - - -def diffusion_distance(time, sigma, mu=mu_0): - """ - Distance at which the signal amplitude is largest for a given time after - shut off. Also referred to as the peak distance - - See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html - - - """ - return np.sqrt(2*time/(mu*sigma)) - - ############################################################################### # # # Base Classes # @@ -113,22 +87,6 @@ def cross_orientation(self, xyz): return np.cross(xyz, orientation) -class BaseTDEM(BaseEM): - - time = properties.Float( - "time after shut-off at which we are evaluating the fields (s)", - - required=True - ) - - 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) - - class BaseElectricDipole(BaseDipole): """ Base class for electric current dipoles diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 768eafe4..f40483e5 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -11,6 +11,11 @@ from .base import BaseElectricDipole, BaseMagneticDipole, BaseEM from .. import spatial +__all__ = [ + 'omega', 'wave_number', 'skin_depth', 'sigma_hat', + 'ElectricDipoleWholeSpace' +] + ############################################################################### # # @@ -263,7 +268,15 @@ def current_density(self, xyz): pass def magnetic_field(self, xyz): - pass + dxyz = self.vector_distance(xyz) + r = self.distance(xyz) + r = spatial.repeat_scalar(r) + kr = self.wave_number*r + + front = self.moment / (4. * np.pi * r*3) * np.exp(-1j * kr) + symmetric_term = ( + + ) def magnetic_flux_density(self, xyz): pass diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 8b0e70ca..8b4055bb 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -3,86 +3,184 @@ from __future__ import print_function from __future__ import unicode_literals -from .base import BaseElectricDipole, BaseTDEM - from scipy.constants import mu_0, pi, epsilon_0 from scipy.special import erf import numpy as np -import warnings +import properties + +from .base import BaseElectricDipole, BaseEM + + +############################################################################### +# # +# Utility Functions # +# # +############################################################################### +def peak_time(z, sigma, mu=mu_0): + """ + Time at which the maximum signal amplitude is observed at a particular + location for a transient plane wave through a homogeneous medium. -class ElectricDipole(BaseTDEM, BaseElectricDipole): + See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html + :param z float: distance from source (m) + :param sigma float: electrical conductivity (S/m) + :param mu float: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m """ - .. todo: write this independent of source orientation with dot products + return (mu * sigma * z**2)/6. + + +def diffusion_distance(time, sigma, mu=mu_0): """ - def __init__(self, **kwargs): - super(ElectricDipole, self).__init__(**kwargs) + Distance at which the signal amplitude is largest for a given time after + shut off. Also referred to as the peak distance + + See: http://em.geosci.xyz/content/maxwell1_fundamentals/plane_waves_in_homogeneous_media/time/analytic_solution.html + """ + return np.sqrt(2*time/(mu*sigma)) + + +def theta(time, sigma, mu=mu_0): + """ + Analog to wavenumber in the frequency domain. See Ward and Hohmann, 1988 + pages 174-175 + """ + + +############################################################################### +# # +# Classes # +# # +############################################################################### + +class BaseTDEM(BaseEM): + + time = properties.Float( + "time after shut-off at which we are evaluating the fields (s)", + required=True, + default=1e-4 + ) + + 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 np.sqrt(self.mu*self.sigma / (4.*self.time)) + return np.sqrt(self.mu*self.sigma) + + +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}} + + """ def electric_field(self, xyz): - dxyz = self.vector_distance(xyz) - x = dxyz[:, 0] - y = dxyz[:, 1] - z = dxyz[:, 2] + """ + Electric field from an electric dipole - root_pi = np.sqrt(np.pi) - r = self.distance(xyz) - r2 = r**2 - r3 = r**3. + .. math:: - theta_r = self.theta * r - e_n_theta_r_2 = np.exp(-theta_r**2) + \mathbf{E} = \frac{1}{\hat{\sigma}} \nabla \nabla \cdot \mathbf{A} + - i \omega \mu \mathbf{A} - erf_thetat_r = erf(theta_r) + """ + dxyz = self.vector_distance(xyz) + r = self.distance(xyz) + r = spatial.repeat_scalar(r) + thetar = self.theta * r + root_pi = np.sqrt(np.pi) - current_term = ( - (self.current * self.length) / - (4. * np.pi * self.sigma * r3) + front = ( + (self.current * self.length) / (4 * np.pi * self.sigma * r**3) ) symmetric_term = ( - - ( - (4./root_pi) * theta_r**3 + (6./root_pi) * theta_r - ) * e_n_theta_r_2 + - 3 * erf_thetat_r + ( + ( + 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 + ) ) - src_orientation_term = ( - - ( - (4./root_pi) * theta_r**3 + (2./root_pi) * theta_r - ) * e_n_theta_r_2 + - erf_thetat_r - ) + oriented_term = ( + ( + 4./root_pi * thetar**3 + 2./root_pi * thetar + ) * np.exp(-thetar) + + erf(thetar) + ) * np.kron(self.orientation, np.ones((dxyz.shape[0], 1))) - if np.all(self.orientation == np.r_[1., 0., 0.]): - ex = current_term * ( - symmetric_term * (x**2/r2) - src_orientation_term - ) - ey = current_term * ( symmetric_term * (x*y)/r2 ) - ez = current_term * ( symmetric_term * (x*z)/r2 ) + return front * (symmetric_term + oriented_term) - elif np.all(self.orientation == np.r_[0., 1., 0.]): - ey = current_term * ( - symmetric_term * (y**2/r2) - src_orientation_term - ) - ez = current_term * ( symmetric_term * (y*z)/r2 ) - ex = current_term * ( symmetric_term * (y*x)/r2 ) + def current_density(self, xyz): + """ + Current density due to a harmonic electric dipole + """ + return self.sigma * self.electric_field(xyz) - elif np.all(self.orientation == np.r_[0., 0., 1.]): - ez = current_term * ( - symmetric_term * (z**2/r2) - src_orientation_term + def magnetic_field(self, xyz): + """ + Magnetic field from an electric dipole + """ + dxyz = self.vector_distance(xyz) + r = self.distance(xyz) + 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) ) - ex = current_term * ( symmetric_term * (z*x)/r2 ) - ey = current_term * ( symmetric_term * (z*y)/r2 ) - else: - raise NotImplementedError + ) + + return - front * self.cross_orientation(xyz) / r + + 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 + + def magnetic_flux_density(self, xyz): + """ + Magnetic flux density from an electric dipole + """ + + return self.mu * self.magnetic_field(xyz) + + 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) - return np.c_[ex, ey, ez] - def current_density(self, xyz): - return self.sigma * self.electric_field(xyz) From a6f2458f3d759c1659b3551653a328f9baa3bf16 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 30 Jun 2017 11:44:50 -0700 Subject: [PATCH 22/41] clean up tdem dipole implementation --- geoana/em/tdem.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 8b4055bb..846ff6bb 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -9,6 +9,7 @@ import properties from .base import BaseElectricDipole, BaseEM +from .. import spatial ############################################################################### @@ -109,7 +110,7 @@ def electric_field(self, xyz): ) symmetric_term = ( - ( + - ( ( 4/root_pi * thetar ** 3 + 6/root_pi * thetar ) * np.exp(-thetar**2) + From 093637e060b1ef4a6c71b2cf8f241e6049069513 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 30 Jun 2017 14:00:47 -0700 Subject: [PATCH 23/41] bug fixes in the tdem electric dipole --- docs/conf.py | 4 +++- geoana/em/tdem.py | 8 ++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 765dab16..e5036754 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -303,7 +303,9 @@ def supress_nonlocal_image_warn(): import sphinx.environment - sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn + sphinx.environment.BuildEnvironment.warn_node = ( + _supress_nonlocal_image_warn + ) def _supress_nonlocal_image_warn(self, msg, node, **kwargs): diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 846ff6bb..427ce5be 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -72,7 +72,7 @@ def diffusion_distance(self): @property def theta(self): - return np.sqrt(self.mu*self.sigma) + return np.sqrt(self.mu*self.sigma/(4.*self.time)) class ElectricDipoleWholeSpace(BaseElectricDipole, BaseTDEM): @@ -110,8 +110,8 @@ def electric_field(self, xyz): ) symmetric_term = ( - - ( - ( + ( + - ( 4/root_pi * thetar ** 3 + 6/root_pi * thetar ) * np.exp(-thetar**2) + 3 * erf(thetar) @@ -123,7 +123,7 @@ def electric_field(self, xyz): oriented_term = ( ( 4./root_pi * thetar**3 + 2./root_pi * thetar - ) * np.exp(-thetar) + + ) * np.exp(-thetar**2) - erf(thetar) ) * np.kron(self.orientation, np.ones((dxyz.shape[0], 1))) From 1fc2c200ee4e8de3aaf918e032b2ef9e79e8204c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 30 Jun 2017 14:13:44 -0700 Subject: [PATCH 24/41] cleanup in the docs, install testing packages from requirements --- .travis.yml | 1 - README.rst | 3 ++- docs/index.rst | 13 +++---------- 3 files changed, 5 insertions(+), 12 deletions(-) diff --git a/.travis.yml b/.travis.yml index 9aff16c8..408f3278 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,7 +12,6 @@ install: - pip install --upgrade pip; - pip install -r requirements.txt - pip install -r requirements_dev.txt - - git clone https://github.com/3ptscience/vectormath; - export PYTHONPATH=$PWD/geoana:$PWD/vectormath:$PYTHONPATH # Run test diff --git a/README.rst b/README.rst index 4698ac22..0fd762a5 100644 --- a/README.rst +++ b/README.rst @@ -13,7 +13,7 @@ geoana .. image:: https://travis-ci.org/simpeg/geoana.svg?branch=master :target: https://travis-ci.org/simpeg/geoana - :alt: travis status + :alt: Travis status .. image:: https://www.quantifiedcode.com/api/v1/project/cb381a23f09245bd855f86eac295d8ec/badge.svg :target: https://www.quantifiedcode.com/app/project/cb381a23f09245bd855f86eac295d8ec @@ -23,5 +23,6 @@ geoana :target: https://www.codacy.com/app/lindseyheagy/geoana?utm_source=github.com&utm_medium=referral&utm_content=simpeg/geoana&utm_campaign=Badge_Grade :alt: codacy status + Interactive geoscience (mostly) analytic functions. diff --git a/docs/index.rst b/docs/index.rst index 50ec6668..e2bf9dfc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,8 +1,9 @@ -.. _index: + .. include:: ../README.rst -Contents: + +**Contents:** .. toctree:: :maxdepth: 2 @@ -10,11 +11,3 @@ Contents: content/earthquake content/em - -Indices and tables -================== - -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` - From 3ab34daa6c0384a4d8eb14cbebdbfac2562fbfac Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 30 Jun 2017 14:33:52 -0700 Subject: [PATCH 25/41] linkcheck tuning --- docs/conf.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/conf.py b/docs/conf.py index e5036754..a457d1c3 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -119,6 +119,13 @@ html_theme = 'sphinx_rtd_theme' html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] + +linkcheck_ignore = [ +] + +linkcheck_retries = 3 +linkcheck_timeout = 500 + # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. From a0c8f66605ccb6f494bd44cfcb14e36d2dd47727 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 1 Jul 2017 10:14:53 -0700 Subject: [PATCH 26/41] use https for targets --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 0fd762a5..16a10983 100644 --- a/README.rst +++ b/README.rst @@ -4,7 +4,7 @@ geoana ====== .. image:: https://readthedocs.org/projects/geoana/badge/?version=latest - :target: http://geoana.readthedocs.io/en/latest/?badge=latest + :target: https://geoana.readthedocs.io/en/latest/?badge=latest :alt: Documentation Status .. image:: https://img.shields.io/github/license/simpeg/geoana.svg From 70cf98a3b2d57e98b85d02ee2aa451352fba2230 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 1 Jul 2017 10:29:45 -0700 Subject: [PATCH 27/41] install using setup.py --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 408f3278..84d7d8b3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,7 +12,7 @@ install: - pip install --upgrade pip; - pip install -r requirements.txt - pip install -r requirements_dev.txt - - export PYTHONPATH=$PWD/geoana:$PWD/vectormath:$PYTHONPATH + - cd geoana; python setup.py install; cd .. # Run test script: From 559bbaba149d1d3ce6f20caf3ffc651f0c2bdf55 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 1 Jul 2017 10:55:24 -0700 Subject: [PATCH 28/41] add rtd badge to linkcheck ignore --- docs/conf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/conf.py b/docs/conf.py index a457d1c3..6b90ffce 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -121,6 +121,7 @@ linkcheck_ignore = [ + 'https://readthedocs.org/projects/geoana/badge/?version=latest' ] linkcheck_retries = 3 From eb29ab01cb4d6390268bd207ff66b272313f57c5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 1 Jul 2017 11:19:15 -0700 Subject: [PATCH 29/41] remove linkcheck --- tests/test_docs.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/test_docs.py b/tests/test_docs.py index fe31514c..01ab9b1e 100644 --- a/tests/test_docs.py +++ b/tests/test_docs.py @@ -43,14 +43,14 @@ def test_html(self): ) assert check == 0 - def test_linkcheck(self): - check = subprocess.call( - ['sphinx-build', '-nW', '-b', 'linkcheck', '-d', - '{}'.format(self.doctrees_dir), - '{}'.format(self.docs_dir), - '{}'.format(self.build_dir)] - ) - assert check == 0 + # def test_linkcheck(self): + # check = subprocess.call( + # ['sphinx-build', '-nW', '-b', 'linkcheck', '-d', + # '{}'.format(self.doctrees_dir), + # '{}'.format(self.docs_dir), + # '{}'.format(self.build_dir)] + # ) + # assert check == 0 if __name__ == '__main__': unittest.main() From 8cede7894d2c908bc2f36a79cfc04fb637857eab Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 3 Jul 2017 08:35:39 -0700 Subject: [PATCH 30/41] add jupyter to conda env --- appveyor.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/appveyor.yml b/appveyor.yml index 9906bc68..a909dec3 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -36,7 +36,7 @@ install: - conda config --set always_yes yes --set changeps1 no - conda update -q conda - conda info -a - - "conda create -q -n test-environment python=%PYTHON% numpy scipy matplotlib cython ipython pillow wheel" + - "conda create -q -n test-environment python=%PYTHON% numpy scipy matplotlib cython jupyter ipython pillow wheel" - activate test-environment - pip install -r requirements_dev.txt - python setup.py install From 2b3b7cd2a7202e1b0d547d9d18ee27248f1c8202 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 3 Jul 2017 08:43:12 -0700 Subject: [PATCH 31/41] add nose to dev requirements --- requirements_dev.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/requirements_dev.txt b/requirements_dev.txt index 23d1fb05..ff204b34 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -2,3 +2,4 @@ discretize sphinx_rtd_theme SimPEG sphinx +nose From 112397598f8c37cd32e8f66e4b96a2fcb7d6a695 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 3 Jul 2017 09:09:12 -0700 Subject: [PATCH 32/41] cleanup unnecessary imports --- geoana/em/fdem.py | 3 ++- geoana/em/tdem.py | 2 +- geoana/spatial.py | 1 - 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index f40483e5..2b06697e 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -248,7 +248,7 @@ def magnetic_field(self, xyz): self.current * self.length / (4 * np.pi * r**2) * (1j * kr + 1) * np.exp(-1j * kr) ) - return - front * self.cross_orientation(xyz) / r + return - front * self.cross_orientation(dxyz) / r def magnetic_flux_density(self, xyz): """ @@ -268,6 +268,7 @@ def current_density(self, xyz): pass def magnetic_field(self, xyz): + pass dxyz = self.vector_distance(xyz) r = self.distance(xyz) r = spatial.repeat_scalar(r) diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 427ce5be..1870c433 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -3,7 +3,7 @@ from __future__ import print_function from __future__ import unicode_literals -from scipy.constants import mu_0, pi, epsilon_0 +from scipy.constants import mu_0, pi from scipy.special import erf import numpy as np import properties diff --git a/geoana/spatial.py b/geoana/spatial.py index 07d429db..04ff604f 100644 --- a/geoana/spatial.py +++ b/geoana/spatial.py @@ -4,7 +4,6 @@ from __future__ import unicode_literals import numpy as np -import properties def vector_distance(xyz, origin=np.r_[0., 0., 0.]): From 96b94cf3e9d13425717b7af0d0f2c3107c104400 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 16 Nov 2017 17:39:50 -0300 Subject: [PATCH 33/41] cleanup in the electric dipole, testing e-field from magnetic dipole --- geoana/__init__.py | 1 + geoana/em/fdem.py | 127 +++++-- geoana/em/tdem.py | 3 +- geoana/utils.py | 105 ++++++ tests/test_em_fdem_magnetic_dipole.py | 524 ++++++++++++++++++++++---- 5 files changed, 643 insertions(+), 117 deletions(-) create mode 100644 geoana/utils.py diff --git a/geoana/__init__.py b/geoana/__init__.py index ae43d999..03a8faf9 100644 --- a/geoana/__init__.py +++ b/geoana/__init__.py @@ -1,5 +1,6 @@ from . import earthquake from . import em +from . import utils __version__ = '0.0.1' __author__ = 'SimPEG developers' diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 2b06697e..73671a74 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -12,8 +12,8 @@ from .. import spatial __all__ = [ - 'omega', 'wave_number', 'skin_depth', 'sigma_hat', - 'ElectricDipoleWholeSpace' + 'omega', 'wavenumber', 'skin_depth', 'sigma_hat', + 'ElectricDipoleWholeSpace', 'MagneticDipoleWholeSpace' ] @@ -36,8 +36,8 @@ def omega(frequency): return 2*np.pi*frequency -def wave_number( - frequency, sigma, mu=mu_0, epsilon=epsilon_0, quasi_static=False +def wavenumber( + frequency, sigma, mu=mu_0, epsilon=epsilon_0, quasistatic=False ): """ Wavenumber of an electromagnetic wave in a medium with constant physical @@ -51,12 +51,12 @@ def wave_number( :param float sigma: electrical conductivity (S/m) :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 quasi_static: use the quasi-static assumption? Default: False + :param bool quasistatic: use the quasi-static assumption? Default: False """ w = omega(frequency) - if quasi_static is True: + if quasistatic is True: return np.sqrt(-1j * w * mu * sigma) - return np.sqrt(w**2. * mu * epsilon - 1j * w * mu * sigma) + return np.sqrt(w**2 * mu * epsilon - 1j * w * mu * sigma) def skin_depth(frequency, sigma, mu=mu_0): @@ -76,7 +76,7 @@ def skin_depth(frequency, sigma, mu=mu_0): return np.sqrt(2./(omega*sigma*mu)) -def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasi_static=False): +def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasistatic=False): """ conductivity with displacement current contribution @@ -87,9 +87,9 @@ def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasi_static=False): :param (float, numpy.array) frequency: frequency (Hz) :param float sigma: electrical conductivity (S/m) :param float epsilon: dielectric permittivity. Default :math:`\varepsilon_0` - :param bool quasi_static: use the quasi-static assumption? Default: False + :param bool quasistatic: use the quasi-static assumption? Default: False """ - if quasi_static is True: + if quasistatic is True: return sigma return sigma + 1j*omega(frequency)*epsilon @@ -106,11 +106,11 @@ class BaseFDEM(BaseEM): """ frequency = properties.Float( "Source frequency (Hz)", - default=1, + default=1., min=0.0 ) - quasi_static = properties.Bool( + quasistatic = properties.Bool( "Use the quasi-static approximation and ignore displacement current?", default=False ) @@ -137,11 +137,12 @@ def sigma_hat(self): """ return sigma_hat( - self.frequency, self.sigma, self.epsilon, self.quasi_static + self.frequency, self.sigma, epsilon=self.epsilon, + quasistatic=self.quasistatic ) @property - def wave_number(self): + def wavenumber(self): """ Wavenumber of an electromagnetic wave in a medium with constant physical properties @@ -150,7 +151,10 @@ def wave_number(self): k = \sqrt{\omega**2 \mu \varepsilon - i \omega \mu \sigma} """ - return wave_number(self.frequency, self.sigma, self.mu) + return wavenumber( + self.frequency, self.sigma, mu=self.mu, epsilon=self.epsilon, + quasistatic=self.quasistatic + ) @property def skin_depth(self): @@ -162,7 +166,7 @@ def skin_depth(self): \sqrt{\\frac{2}{\omega \sigma \mu}} """ - return skin_depth(self.frequency, self.sigma, self.mu) + return skin_depth(self.frequency, self.sigma, mu=self.mu) class ElectricDipoleWholeSpace(BaseElectricDipole, BaseFDEM): @@ -190,7 +194,7 @@ def vector_potential(self, xyz): r = self.distance(xyz) a = ( (self.current * self.length) / (4*np.pi*r) * - np.exp(-i*self.wave_number*r) + np.exp(-i*self.wavenumber*r) ) a = np.kron(np.ones(1, 3), np.atleast_2d(a).T) return self.dot_orientation(a) @@ -206,23 +210,23 @@ def electric_field(self, xyz): """ dxyz = self.vector_distance(xyz) - r = self.distance(xyz) - r = spatial.repeat_scalar(r) - kr = self.wave_number * r + r = spatial.repeat_scalar(self.distance(xyz)) + kr = self.wavenumber * r + ikr = 1j * kr - front = ( + front_term = ( (self.current * self.length) / (4 * np.pi * self.sigma * r**3) * - np.exp(-1j * kr) + np.exp(-ikr) ) symmetric_term = ( spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz * - (-kr**2 + 3*1j*kr + 3) / r**2 + (-kr**2 + 3*ikr+ 3) / r**2 ) oriented_term = ( - (kr**2 - 1j*kr - 1) * + (kr**2 - ikr - 1) * np.kron(self.orientation, np.ones((dxyz.shape[0], 1))) ) - return front * (symmetric_term + oriented_term) + return front_term * (symmetric_term + oriented_term) def current_density(self, xyz): """ @@ -240,15 +244,15 @@ def magnetic_field(self, xyz): """ dxyz = self.vector_distance(xyz) - r = self.distance(xyz) - r = spatial.repeat_scalar(r) - kr = self.wave_number * r + r = spatial.repeat_scalar(self.distance(xyz)) + kr = self.wavenumber * r + ikr = 1j*kr - front = ( - self.current * self.length / (4 * np.pi * r**2) * (1j * kr + 1) * - np.exp(-1j * kr) + front_term = ( + self.current * self.length / (4 * np.pi * r**2) * (ikr + 1) * + np.exp(-ikr) ) - return - front * self.cross_orientation(dxyz) / r + return -front_term * self.cross_orientation(dxyz) / r def magnetic_flux_density(self, xyz): """ @@ -261,23 +265,68 @@ class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseFDEM): """ Harmonic magnetic dipole in a whole space. """ + + 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) + def electric_field(self, xyz): - pass + """ + 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 + + front_term = ( + (1j * self.omega * self.mu * self.moment) / (4. * np.pi * r**2) * + (1j * kr + 1) * np.exp(-1j * kr) + ) + return front_term * self.cross_orientation(dxyz) / r def current_density(self, xyz): - pass + """ + Current density from a magnetic dipole in a wholespace + """ + return self.sigma * self.electric_field(xyz) def magnetic_field(self, xyz): - pass + """ + Magnetic field due to a magnetic dipole in a wholespace + """ dxyz = self.vector_distance(xyz) r = self.distance(xyz) r = spatial.repeat_scalar(r) - kr = self.wave_number*r + kr = self.wavenumber*r + ikr = 1j*kr - front = self.moment / (4. * np.pi * r*3) * np.exp(-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) + def magnetic_flux_density(self, xyz): - pass + """ + Magnetic flux density due to a magnetic dipole in a wholespace + """ + return self.mu * self.magnetic_field(xyz) diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index 1870c433..a5c78460 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -47,6 +47,7 @@ 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)) ############################################################################### @@ -72,7 +73,7 @@ def diffusion_distance(self): @property def theta(self): - return np.sqrt(self.mu*self.sigma/(4.*self.time)) + return theta(self.time, self.sigma, self.mu) class ElectricDipoleWholeSpace(BaseElectricDipole, BaseTDEM): diff --git a/geoana/utils.py b/geoana/utils.py new file mode 100644 index 00000000..f9e69c95 --- /dev/null +++ b/geoana/utils.py @@ -0,0 +1,105 @@ +import numpy as np + + +def mkvc(x, numDims=1): + """Creates a vector with the number of dimension specified + + e.g.:: + + a = np.array([1, 2, 3]) + + mkvc(a, 1).shape + > (3, ) + + mkvc(a, 2).shape + > (3, 1) + + mkvc(a, 3).shape + > (3, 1, 1) + + """ + if type(x) == np.matrix: + x = np.array(x) + + if hasattr(x, 'tovec'): + x = x.tovec() + + assert isinstance(x, np.ndarray), "Vector must be a numpy array" + + if numDims == 1: + return x.flatten(order='F') + elif numDims == 2: + return x.flatten(order='F')[:, np.newaxis] + elif numDims == 3: + return x.flatten(order='F')[:, np.newaxis, np.newaxis] + + +def ndgrid(*args, **kwargs): + """ + Form tensorial grid for 1, 2, or 3 dimensions. + + Returns as column vectors by default. + + To return as matrix input: + + ndgrid(..., vector=False) + + The inputs can be a list or separate arguments. + + e.g.:: + + a = np.array([1, 2, 3]) + b = np.array([1, 2]) + + XY = ndgrid(a, b) + > [[1 1] + [2 1] + [3 1] + [1 2] + [2 2] + [3 2]] + + X, Y = ndgrid(a, b, vector=False) + > X = [[1 1] + [2 2] + [3 3]] + > Y = [[1 2] + [1 2] + [1 2]] + + """ + + # Read the keyword arguments, and only accept a vector=True/False + vector = kwargs.pop('vector', True) + assert type(vector) == bool, "'vector' keyword must be a bool" + assert len(kwargs) == 0, "Only 'vector' keyword accepted" + + # you can either pass a list [x1, x2, x3] or each seperately + if type(args[0]) == list: + xin = args[0] + else: + xin = args + + # Each vector needs to be a numpy array + assert np.all( + [isinstance(x, np.ndarray) for x in xin] + ), "All vectors must be numpy arrays." + + if len(xin) == 1: + return xin[0] + elif len(xin) == 2: + XY = np.broadcast_arrays(mkvc(xin[1], 1), mkvc(xin[0], 2)) + if vector: + X2, X1 = [mkvc(x) for x in XY] + return np.c_[X1, X2] + else: + return XY[1], XY[0] + elif len(xin) == 3: + XYZ = np.broadcast_arrays( + mkvc(xin[2], 1), mkvc(xin[1], 2), mkvc(xin[0], 3) + ) + if vector: + X3, X2, X1 = [mkvc(x) for x in XYZ] + return np.c_[X1, X2, X3] + else: + return XYZ[2], XYZ[1], XYZ[0] diff --git a/tests/test_em_fdem_magnetic_dipole.py b/tests/test_em_fdem_magnetic_dipole.py index c7a48181..dc9adf4d 100644 --- a/tests/test_em_fdem_magnetic_dipole.py +++ b/tests/test_em_fdem_magnetic_dipole.py @@ -5,108 +5,478 @@ import unittest import numpy as np +import discretize +from SimPEG import Maps +from SimPEG.EM import FDEM from scipy.constants import mu_0, epsilon_0 from geoana.em import fdem from discretize.utils import ndgrid, asArray_N_x_Dim -def MagneticDipoleFields(srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0): +def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=1., epsr=1., t=0.): + + """ + Computing magnetic fields from Magnetic Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0 * (1+kappa) + epsilon = epsilon_0 * epsr + m = current * loopArea + + omega = lambda f: 2*np.pi*f + + XYZ = discretize.utils.asArray_N_x_Dim(XYZ, 3) + # Check + + dx = XYZ[:, 0]-srcLoc[0] + dy = XYZ[:, 1]-srcLoc[1] + dz = XYZ[:, 2]-srcLoc[2] + + r = np.sqrt(dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt(omega(f)**2. * mu*epsilon - 1j*omega(f)*mu*sig) + + front = m / (4.*np.pi*(r)**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Hx = front*((dx**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) + Hy = front*(dx*dy / r**2)*mid + Hz = front*(dx*dz / r**2)*mid + return Hx, Hy, Hz + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Hy = front * ((dy**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) + Hz = front * (dy*dz / r**2)*mid + Hx = front * (dy*dx / r**2)*mid + return Hx, Hy, Hz + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Hz = front*((dz**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) + Hx = front*(dz*dx / r**2)*mid + Hy = front*(dz*dy / r**2)*mid + return Hx, Hy, Hz + + +def B_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=1., epsr=1., t=0.): + """ - Calculate the vector potential of a set of magnetic dipoles - at given locations 'ref. ' - - .. math:: - - B = \frac{\mu_0}{4 \pi r^3} \left( \frac{3 \vec{r} (\vec{m} \cdot - \vec{r})}{r^2}) - - \vec{m} - \right) \cdot{\hat{rx}} - - :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) - :param numpy.ndarray obsLoc: Where the potentials will be calculated - (x, y, z) - :param str component: The component to calculate - 'x', 'y', or 'z' - :param numpy.ndarray moment: The vector dipole moment (vertical) - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location + Computing magnetic flux densites from Magnetic Dipole in a Wholespace + TODO: + Add description of parameters """ + mu = mu_0 * (1+kappa) + + Hx, Hy, Hz = H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, loopArea=loopArea, orientation=orientation, kappa=kappa, epsr=epsr) + Bx = mu * Hx + By = mu * Hy + Bz = mu * Hz + return Bx, By, Bz + + +def E_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0., epsr=1., t=0.): + + """ + Computing analytic electric fields from Magnetic Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0 * (1+kappa) + epsilon = epsilon_0 * epsr + m = current * loopArea + + omega = lambda f: 2 * np.pi * f + + XYZ = discretize.utils.asArray_N_x_Dim(XYZ, 3) + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + + front = ((1j * omega(f) * mu * m) / (4.* np.pi * r**2)) * (1j * k * r + 1) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Ey = front * (dz / r) + Ez = front * (-dy / r) + Ex = np.zeros_like(Ey) + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + Ex = front * (-dz / r) + Ez = front * (dx / r) + Ey = np.zeros_like(Ex) + return Ex, Ey, Ez + + elif orientation.upper() == 'Z': + Ex = front * (dy / r) + Ey = front * (-dx / r) + Ez = np.zeros_like(Ex) + return Ex, Ey, Ez - if isinstance(orientation, str): - assert orientation.upper() in ['X', 'Y', 'Z'], ( - "orientation must be 'x', 'y', or 'z' or a vector not {}" - .format(orientation) + +class TestFDEMdipole(unittest.TestCase): + + def test_defaults(self): + TOL = 1e-15 + mdws = fdem.MagneticDipoleWholeSpace() + assert(mdws.sigma == 1.) + assert(mdws.mu == mu_0) + assert(mdws.epsilon == epsilon_0) + assert(np.all(mdws.orientation == np.r_[1., 0., 0.])) + assert(mdws.moment == 1.) + assert(np.all(mdws.location == np.r_[0., 0., 0.])) + assert(mdws.frequency == 1.) + assert(mdws.omega == 2.*np.pi*1.) + assert(mdws.quasistatic is False) + assert np.linalg.norm( + mdws.wavenumber - np.sqrt( + mu_0 * epsilon_0 * (2*np.pi)**2 - 1j * mu_0 * 1. * 2*np.pi + ) + ) <= TOL + assert np.linalg.norm( + mdws.wavenumber**2 - ( + mu_0 * epsilon_0 * (2*np.pi)**2 - 1j * mu_0 * 1. * 2*np.pi + ) + ) <= TOL + + def compare_fields(name, field, ftest): + + def check_component(name, f, ftest): + geoana_norm = np.linalg.norm(f) + test_norm = np.linalg.norm(ftest) + diff = np.linalg.norm(f-ftest) + passed = np.allclose(f, ftest) + print( + "Testing {} ... geoana: {:1.4e}, compare: {:1.4e}, " + "diff: {:1.4e}, passed?: {}".format( + name, geoana_norm, test_norm, diff, passed + ) + ) + return passed + + passed = [] + for i, orientation in enumerate(['x', 'y', 'z']): + for component in ['real', 'imag']: + passed.append(check_component( + orientation + '_' + component, + getattr(field[:, i], component), + getattr(ftest[:, i], component) + )) + return all(passed) + + def magnetic_dipole_b(self, orientation): + sigma = 1e-3 + frequency = 1. + mdws = fdem.MagneticDipoleWholeSpace( + orientation=orientation, + sigma=sigma, + frequency=frequency + ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + xyz = discretize.utils.ndgrid([x, y, z]) + + # srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0 + + # btest = [MagneticDipoleFields( + # mdws.location, xyz, rx_orientation, + # orientation=orientation.upper() + # ) for rx_orientation in ["x", "y", "z"]] + + bxtest, bytest, bztest = B_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency ) - elif (not np.allclose(np.r_[1., 0., 0.], orientation) or - not np.allclose(np.r_[0., 1., 0.], orientation) or - not np.allclose(np.r_[0., 0., 1.], orientation)): - warnings.warn( - 'Arbitrary trasnmitter orientations ({}) not thouroughly tested ' - 'Pull request on a test anyone? bueller?'.format(orientation) + + b = mdws.magnetic_flux_density(xyz) + print( + "\n\nTesting Magnetic Dipole B: {} orientation\n".format(orientation) + ) + + passed = self.compare_fields(b, np.vstack([bxtest, bytest, bztest]).T) + self.assertTrue(passed) + + def magnetic_dipole_e(self, orientation): + sigma = 1e-2 + frequency = 1 + mdws = fdem.MagneticDipoleWholeSpace( + orientation=orientation, + sigma=sigma, + frequency=frequency ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + xyz = discretize.utils.ndgrid([x, y, z]) + + extest, eytest, eztest = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation=orientation + ) + + e = mdws.electric_field(xyz) + print( + "\n\nTesting Magnetic Dipole E: {} orientation\n".format(orientation) + ) + + passed = self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + self.assertTrue(passed) + + def test_magnetic_dipole_x_b(self): + self.magnetic_dipole_b("x") + + def test_magnetic_dipole_y_b(self): + self.magnetic_dipole_b("y") + + def test_magnetic_dipole_z_b(self): + self.magnetic_dipole_b("z") + + def test_magnetic_dipole_tilted_b(self): + + orientation = np.random.rand(3) + orientation = orientation / np.linalg.norm(orientation) + + mdws = fdem.MagneticDipoleWholeSpace( + orientation=orientation + ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + + xyz = discretize.utils.ndgrid([x, y, z]) - if isinstance(component, str): - assert component.upper() in ['X', 'Y', 'Z'], ( - "component must be 'x', 'y', or 'z' or a vector not {}" - .format(component) + bxtest0, bytest0, bztest0 = B_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='X' ) - elif (not np.allclose(np.r_[1., 0., 0.], component) or - not np.allclose(np.r_[0., 1., 0.], component) or - not np.allclose(np.r_[0., 0., 1.], component)): - warnings.warn( - 'Arbitrary receiver orientations ({}) not thouroughly tested ' - 'Pull request on a test anyone? bueller?' - ).format(component) + bxtest1, bytest1, bztest1 = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='Y' + ) + bxtest2, bytest2, bztest2 = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='Z' + ) + + bxtest = ( + orientation[0]*bxtest0 + orientation[1]*bxtest1 + orientation[2]*bxtest2 + ) + bytest = ( + orientation[0]*bytest0 + orientation[1]*bytest1 + orientation[2]*bytest2 + ) + bztest = ( + orientation[0]*bztest0 + orientation[1]*bztest1 + orientation[2]*bztest2 + ) + + b = mdws.magnetic_flux_density(xyz) + print( + "\n\nTesting Magnetic Dipole B: {} orientation\n".format("45 degree") + ) + + self.compare_fields(b, np.vstack([bxtest, bytest, bztest]).T) + + def test_magnetic_dipole_x_e(self): + self.magnetic_dipole_e("x") + + def test_magnetic_dipole_y_e(self): + self.magnetic_dipole_e("y") + + def test_magnetic_dipole_z_e(self): + self.magnetic_dipole_e("z") + + def test_magnetic_dipole_tilted_e(self): + + orientation = np.random.rand(3) + orientation = orientation / np.linalg.norm(orientation) + + mdws = fdem.MagneticDipoleWholeSpace( + orientation=orientation + ) + x = np.linspace(-20., 20., 50) + y = np.linspace(-30., 30., 50) + z = np.linspace(-40., 40., 50) + + xyz = discretize.utils.ndgrid([x, y, z]) + + extest0, eytest0, eztest0 = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='X' + ) + extest1, eytest1, eztest1 = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='Y' + ) + extest2, eytest2, eztest2 = E_from_MagneticDipoleWholeSpace( + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation='Z' + ) + + extest = ( + orientation[0]*extest0 + orientation[1]*extest1 + orientation[2]*extest2 + ) + eytest = ( + orientation[0]*eytest0 + orientation[1]*eytest1 + orientation[2]*eytest2 + ) + eztest = ( + orientation[0]*eztest0 + orientation[1]*eztest1 + orientation[2]*eztest2 + ) + + e = mdws.electric_field(xyz) + print( + "\n\nTesting Magnetic Dipole E: {} orientation\n".format("45 degree") + ) + + self.compare_fields(e, np.vstack([extest, eytest, eztest]).T) + + +class TestFDEMdipole_SimPEG(unittest.TestCase): + + tol = 1e-1 # error must be an order of magnitude smaller than results + + # put the source at the center + + def getFaceSrc(self, mesh): + csx = mesh.hx.min() + csz = mesh.hz.min() + srcInd = ( + (mesh.gridFz[:, 0] < csx) & + (mesh.gridFz[:, 2] < csz/2.) & (mesh.gridFz[:, 2] > -csz/2.) + ) + + src_vecz = np.zeros(mesh.nFz, dtype=complex) + src_vecz[srcInd] = 1. + + return np.hstack( + [np.zeros(mesh.vnF[:2].sum(), dtype=complex), src_vecz] + ) + + def getProjections(self, mesh): + ignore_inside_radius = 5*mesh.hx.min() + ignore_outside_radius = 40*mesh.hx.min() + + def ignoredGridLocs(grid): + return ( + ( + grid[:, 0]**2 + grid[:, 1]**2 + grid[:, 2]**2 < + ignore_inside_radius**2 + ) | ( + grid[:, 0]**2 + grid[:, 1]**2 + grid[:, 2]**2 > + ignore_outside_radius**2 + ) + ) + + # Faces + ignore_me_Fx = ignoredGridLocs(mesh.gridFx) + ignore_me_Fz = ignoredGridLocs(mesh.gridFz) + ignore_me = np.hstack([ignore_me_Fx, ignore_me_Fz]) + keep_me = np.array(~ignore_me, dtype=float) + Pf = discretize.utils.sdiag(keep_me) + + # Edges + ignore_me_Ey = ignoredGridLocs(mesh.gridEy) + keep_me_Ey = np.array(~ignore_me_Ey, dtype=float) + Pe = discretize.utils.sdiag(keep_me_Ey) + + return Pf, Pe + + def test_b_dipole_v_SimPEG(self): + + def compare_w_SimPEG(name, geoana, simpeg): + + norm_geoana = np.linalg.norm(geoana) + norm_simpeg = np.linalg.norm(simpeg) + diff = np.linalg.norm(geoana - simpeg) + passed = diff < self.tol * 0.5 * (norm_geoana + norm_simpeg) + print( + " {} ... geoana: {:1.4e}, SimPEG: {:1.4e}, diff: {:1.4e}, " + "passed?: {}".format( + name, norm_geoana, norm_simpeg, diff, passed + ) + ) + + return passed + + print("\n\nComparing Magnetic dipole with SimPEG") + + sigma_back = 1e-1 + freqs = np.r_[10., 100.] + + csx, ncx, npadx = 0.5, 50, 30 + ncy = 1 + csz, ncz, npadz = 0.5, 50, 30 + + hx = discretize.utils.meshTensor( + [(csx, ncx), (csx, npadx, 1.2)] + ) + hy = 2*np.pi / ncy * np.ones(ncy) + hz = discretize.utils.meshTensor( + [(csz, npadz, -1.2), (csz, ncz), (csz, npadz, 1.2)] + ) + + mesh = discretize.CylMesh([hx, hy, hz], x0='00C') + + s_m = self.getFaceSrc(mesh) + prob = FDEM.Problem3D_b(mesh, sigmaMap=Maps.IdentityMap(mesh)) + srcList = [FDEM.Src.RawVec_m([], f, s_m) for f in freqs] + survey = FDEM.Survey(srcList) - if isinstance(orientation, str): - orientation = orientationDict[orientation.upper()] + prob.pair(survey) - if isinstance(component, str): - component = orientationDict[component.upper()] + fields = prob.fields(sigma_back*np.ones(mesh.nC)) - assert np.linalg.norm(orientation, 2) == 1., ( - "orientation must be a unit vector. " - "Use 'moment=X to scale source fields" - ) + moment = np.pi * csx**2 # mesh.hz.min() * - if np.linalg.norm(component, 2) != 1.: - warnings.warn( - 'The magnitude of the receiver component vector is > 1, ' - ' it is {}. The receiver fields will be scaled.' - .format(np.linalg.norm(component, 2)) + mdws = fdem.MagneticDipoleWholeSpace( + sigma=sigma_back, moment=moment, orientation="z" ) - srcLoc = np.atleast_2d(srcLoc) - component = np.atleast_2d(component) - obsLoc = np.atleast_2d(obsLoc) - orientation = np.atleast_2d(orientation) + Pf, Pe = self.getProjections(mesh) - nObs = obsLoc.shape[0] - nSrc = int(srcLoc.size / 3.) + e_passed = [] + b_passed = [] + for i, f in enumerate(freqs): + mdws.frequency = f - # use outer product to construct an array of [x_src, y_src, z_src] + b_xz = [] + for b, component in zip([0, 2], ['x', 'z']): + grid = getattr(mesh, "gridF{}".format(component)) + b_xz.append( + mdws.magnetic_flux_density(grid)[:, b] + ) + b_geoana = np.hstack(b_xz) + e_geoana = mdws.electric_field(mesh.gridEy)[:, 1] - m = moment*orientation.repeat(nObs, axis=0) - B = [] + P_e_geoana = Pe*e_geoana + P_e_simpeg = Pe*discretize.utils.mkvc(fields[srcList[i], 'e']) - for i in range(nSrc): - srcLoc = srcLoc[i, np.newaxis].repeat(nObs, axis=0) - rx = component.repeat(nObs, axis=0) - dR = obsLoc - srcLoc - r = np.sqrt((dR**2).sum(axis=1)) + P_b_geoana = Pf*b_geoana + P_b_simpeg = Pf*discretize.utils.mkvc(fields[srcList[i], 'b']) - # mult each element and sum along the axis (vector dot product) - m_dot_dR_div_r2 = (m * dR).sum(axis=1) / (r**2) + print("Testing {} Hz".format(f)) - # multiply the scalar m_dot_dR by the 3D vector r - rvec_m_dot_dR_div_r2 = np.vstack([m_dot_dR_div_r2 * dR[:, i] for - i in range(3)]).T - inside = (3. * rvec_m_dot_dR_div_r2) - m + for comp in ['real', 'imag']: + e_passed.append(compare_w_SimPEG( + 'E {}'.format(comp), + getattr(P_e_geoana, comp), + getattr(P_e_simpeg, comp) + )) + b_passed.append(compare_w_SimPEG( + 'B {}'.format(comp), + getattr(P_b_geoana, comp), + getattr(P_b_simpeg, comp) + )) + assert(all(e_passed)) + assert(all(b_passed)) - # dot product with rx orientation - inside_dot_rx = (inside * rx).sum(axis=1) - front = (mu/(4.* pi * r**3)) - B.append(Utils.mkvc(front * inside_dot_rx)) +if __name__ == '__main__': + unittest.main() - return np.vstack(B).T From 21fa67b554a76ab1143d44b5bb5106cc05fd7f19 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 28 Nov 2017 14:40:40 -0200 Subject: [PATCH 34/41] clean up mag dipole code, add tests against independently developed analytic, compare with SimPEG --- geoana/em/fdem.py | 8 +-- tests/test_em_fdem_magnetic_dipole.py | 84 +++++++++++++++------------ 2 files changed, 50 insertions(+), 42 deletions(-) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index 73671a74..ae958cf5 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -220,7 +220,7 @@ def electric_field(self, xyz): ) symmetric_term = ( spatial.repeat_scalar(self.dot_orientation(dxyz)) * dxyz * - (-kr**2 + 3*ikr+ 3) / r**2 + (-kr**2 + 3*ikr + 3) / r**2 ) oriented_term = ( (kr**2 - ikr - 1) * @@ -290,10 +290,11 @@ def electric_field(self, xyz): 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) * - (1j * kr + 1) * np.exp(-1j * kr) + (ikr + 1) * np.exp(-ikr) ) return front_term * self.cross_orientation(dxyz) / r @@ -308,8 +309,7 @@ def magnetic_field(self, xyz): Magnetic field due to a magnetic dipole in a wholespace """ dxyz = self.vector_distance(xyz) - r = self.distance(xyz) - r = spatial.repeat_scalar(r) + r = spatial.repeat_scalar(self.distance(xyz)) kr = self.wavenumber*r ikr = 1j*kr diff --git a/tests/test_em_fdem_magnetic_dipole.py b/tests/test_em_fdem_magnetic_dipole.py index dc9adf4d..79e65e89 100644 --- a/tests/test_em_fdem_magnetic_dipole.py +++ b/tests/test_em_fdem_magnetic_dipole.py @@ -14,17 +14,25 @@ from discretize.utils import ndgrid, asArray_N_x_Dim -def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=1., epsr=1., t=0.): +def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, epsr=1., t=0.): """ Computing magnetic fields from Magnetic Dipole in a Wholespace TODO: Add description of parameters """ + assert current == 1 + assert loopArea == 1 + assert np.all(srcLoc == np.r_[0., 0., 0.]) + assert kappa == 0 mu = mu_0 * (1+kappa) epsilon = epsilon_0 * epsr m = current * loopArea + assert m == 1 + assert mu == mu_0 + assert epsilon == epsilon_0 + omega = lambda f: 2*np.pi*f XYZ = discretize.utils.asArray_N_x_Dim(XYZ, 3) @@ -36,33 +44,32 @@ def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1. r = np.sqrt(dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt(omega(f)**2. * mu*epsilon - 1j*omega(f)*mu*sig) + k = np.sqrt(omega(f)**2. * mu * epsilon - 1j*omega(f)*mu*sig) - front = m / (4.*np.pi*(r)**3) * np.exp(-1j*k*r) - mid = -k**2 * r**2 + 3*1j*k*r + 3 + front = m / (4. * np.pi * r**3) * np.exp(-1j*k*r) + mid = - k**2 * r**2 + 3*1j*k*r + 3 if orientation.upper() == 'X': Hx = front*((dx**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) Hy = front*(dx*dy / r**2)*mid Hz = front*(dx*dz / r**2)*mid - return Hx, Hy, Hz elif orientation.upper() == 'Y': # x--> y, y--> z, z-->x Hy = front * ((dy**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) Hz = front * (dy*dz / r**2)*mid Hx = front * (dy*dx / r**2)*mid - return Hx, Hy, Hz elif orientation.upper() == 'Z': # x --> z, y --> x, z --> y Hz = front*((dz**2 / r**2)*mid + (k**2 * r**2 - 1j*k*r - 1.)) Hx = front*(dz*dx / r**2)*mid Hy = front*(dz*dy / r**2)*mid - return Hx, Hy, Hz + + return Hx, Hy, Hz -def B_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=1., epsr=1., t=0.): +def B_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, epsr=1., t=0.): """ Computing magnetic flux densites from Magnetic Dipole in a Wholespace @@ -93,9 +100,9 @@ def E_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1. XYZ = discretize.utils.asArray_N_x_Dim(XYZ, 3) - dx = XYZ[:,0]-srcLoc[0] - dy = XYZ[:,1]-srcLoc[1] - dz = XYZ[:,2]-srcLoc[2] + dx = XYZ[:, 0]-srcLoc[0] + dy = XYZ[:, 1]-srcLoc[1] + dz = XYZ[:, 2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) @@ -173,7 +180,7 @@ def check_component(name, f, ftest): return all(passed) def magnetic_dipole_b(self, orientation): - sigma = 1e-3 + sigma = 1 frequency = 1. mdws = fdem.MagneticDipoleWholeSpace( orientation=orientation, @@ -193,7 +200,8 @@ def magnetic_dipole_b(self, orientation): # ) for rx_orientation in ["x", "y", "z"]] bxtest, bytest, bztest = B_from_MagneticDipoleWholeSpace( - xyz, mdws.location, mdws.sigma, mdws.frequency + xyz, mdws.location, mdws.sigma, mdws.frequency, + orientation=orientation ) b = mdws.magnetic_flux_density(xyz) @@ -257,11 +265,11 @@ def test_magnetic_dipole_tilted_b(self): xyz, mdws.location, mdws.sigma, mdws.frequency, orientation='X' ) - bxtest1, bytest1, bztest1 = E_from_MagneticDipoleWholeSpace( + bxtest1, bytest1, bztest1 = B_from_MagneticDipoleWholeSpace( xyz, mdws.location, mdws.sigma, mdws.frequency, orientation='Y' ) - bxtest2, bytest2, bztest2 = E_from_MagneticDipoleWholeSpace( + bxtest2, bytest2, bztest2 = B_from_MagneticDipoleWholeSpace( xyz, mdws.location, mdws.sigma, mdws.frequency, orientation='Z' ) @@ -343,23 +351,25 @@ class TestFDEMdipole_SimPEG(unittest.TestCase): # put the source at the center - def getFaceSrc(self, mesh): - csx = mesh.hx.min() - csz = mesh.hz.min() - srcInd = ( - (mesh.gridFz[:, 0] < csx) & - (mesh.gridFz[:, 2] < csz/2.) & (mesh.gridFz[:, 2] > -csz/2.) - ) + # def getFaceSrc(self, mesh): + # csx = mesh.hx.min() + # csz = mesh.hz.min() + # srcInd = ( + # (mesh.gridFz[:, 0] < csx) & + # (mesh.gridFz[:, 2] < csz/2.) & (mesh.gridFz[:, 2] > -csz/2.) + # ) - src_vecz = np.zeros(mesh.nFz, dtype=complex) - src_vecz[srcInd] = 1. + # assert srcInd. - return np.hstack( - [np.zeros(mesh.vnF[:2].sum(), dtype=complex), src_vecz] - ) + # src_vecz = np.zeros(mesh.nFz, dtype=complex) + # src_vecz[srcInd] = 1. + + # return np.hstack( + # [np.zeros(mesh.vnF[:2].sum(), dtype=complex), src_vecz] + # ) def getProjections(self, mesh): - ignore_inside_radius = 5*mesh.hx.min() + ignore_inside_radius = 10*mesh.hx.min() ignore_outside_radius = 40*mesh.hx.min() def ignoredGridLocs(grid): @@ -406,34 +416,32 @@ def compare_w_SimPEG(name, geoana, simpeg): print("\n\nComparing Magnetic dipole with SimPEG") - sigma_back = 1e-1 + sigma_back = 1. freqs = np.r_[10., 100.] - csx, ncx, npadx = 0.5, 50, 30 + csx, ncx, npadx = 1, 50, 50 ncy = 1 - csz, ncz, npadz = 0.5, 50, 30 + csz, ncz, npadz = 1, 50, 50 hx = discretize.utils.meshTensor( - [(csx, ncx), (csx, npadx, 1.2)] + [(csx, ncx), (csx, npadx, 1.3)] ) hy = 2*np.pi / ncy * np.ones(ncy) hz = discretize.utils.meshTensor( - [(csz, npadz, -1.2), (csz, ncz), (csz, npadz, 1.2)] + [(csz, npadz, -1.3), (csz, ncz), (csz, npadz, 1.3)] ) mesh = discretize.CylMesh([hx, hy, hz], x0='00C') - s_m = self.getFaceSrc(mesh) - prob = FDEM.Problem3D_b(mesh, sigmaMap=Maps.IdentityMap(mesh)) - srcList = [FDEM.Src.RawVec_m([], f, s_m) for f in freqs] + prob = FDEM.Problem3D_e(mesh, sigmaMap=Maps.IdentityMap(mesh)) + srcList = [FDEM.Src.MagDipole([], loc=np.r_[0., 0., 0.], freq=f) for f in freqs] survey = FDEM.Survey(srcList) prob.pair(survey) fields = prob.fields(sigma_back*np.ones(mesh.nC)) - moment = np.pi * csx**2 # mesh.hz.min() * - + moment = 1. mdws = fdem.MagneticDipoleWholeSpace( sigma=sigma_back, moment=moment, orientation="z" ) From e3ca77fef7d465628d6648758f24d98e8d08889b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 28 Nov 2017 17:47:43 -0200 Subject: [PATCH 35/41] add a vector magnitude --- geoana/spatial.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/geoana/spatial.py b/geoana/spatial.py index 04ff604f..f90e0b03 100644 --- a/geoana/spatial.py +++ b/geoana/spatial.py @@ -6,6 +6,22 @@ import numpy as np +def vector_magnitude(v): + """ + Amplitude of a vector, v. + + :param numpy.array v: vector array :code:`np.r_[x, y, z]`, with shape (n, 3) + """ + + assert (v.shape[1] == 3), ( + "the vector, v, should be npoints by 3. The shape provided is {}".format( + v.shape + ) + ) + + return np.sqrt((v**2).sum(axis=1)) + + def vector_distance(xyz, origin=np.r_[0., 0., 0.]): """ Vector distance of a grid, xyz from an origin origin. @@ -32,7 +48,7 @@ def distance(xyz, origin=np.r_[0., 0., 0.]): :param numpy.array origin: origin (default: [0., 0., 0.]) """ dxyz = vector_distance(xyz, origin) - return np.sqrt((dxyz**2).sum(axis=1)) + return vector_magnitude(dxyz) def vector_dot(xyz, vector): From f44b24cc0016c7f3b823a55a15a6fae123dbb48d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 28 Nov 2017 17:48:17 -0200 Subject: [PATCH 36/41] add an example plotting the magnetic flux density and electric field due to an electric dipole in a wholespace --- examples/README.txt | 2 + examples/em/README.txt | 2 + examples/em/fdem_ElectricDipoleWholeSpace.py | 113 +++++++++++++++++++ 3 files changed, 117 insertions(+) create mode 100644 examples/README.txt create mode 100644 examples/em/README.txt create mode 100644 examples/em/fdem_ElectricDipoleWholeSpace.py diff --git a/examples/README.txt b/examples/README.txt new file mode 100644 index 00000000..9f935a38 --- /dev/null +++ b/examples/README.txt @@ -0,0 +1,2 @@ +Examples +******** diff --git a/examples/em/README.txt b/examples/em/README.txt new file mode 100644 index 00000000..8ccc45e1 --- /dev/null +++ b/examples/em/README.txt @@ -0,0 +1,2 @@ +Electromagnetics +**************** diff --git a/examples/em/fdem_ElectricDipoleWholeSpace.py b/examples/em/fdem_ElectricDipoleWholeSpace.py new file mode 100644 index 00000000..ccd0dd2a --- /dev/null +++ b/examples/em/fdem_ElectricDipoleWholeSpace.py @@ -0,0 +1,113 @@ +""" +Electric Dipole in a Whole Space: Frequency Domain +================================================== + +In this example, we plot electric and magnetic flux density due to an electric +dipole in a whole space. Note that you can also examine the current density +and magnetic field. + +We can vary the conductivity, magnetic permeability and dielectric permittivity +of the wholespace, the frequency of the source and whether or not the +quasistatic assumption is imposed. + +""" + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +from scipy.constants import mu_0, epsilon_0 + +from geoana import utils, spatial +from geoana.em import fdem + +# define frequencies that we want to look at +frequencies = np.logspace(0, 4, 3) + +# Build the electric dipole object +edipole = fdem.ElectricDipoleWholeSpace( + sigma=1., # conductivity of 1 S/m + mu=mu_0, # permeability of free space (this is the default) + epsilon=epsilon_0, # permittivity of free space (this is the default) + location=np.r_[0., 0., 0.], # location of the dipole + orientation='Z', # vertical dipole (can also be a unit-vector) + quasistatic=False # don't use the quasistatic assumption +) + +# construct a grid where we want to plot electric fields +x = np.linspace(-50, 50, 100) +z = np.linspace(-50, 50, 100) +xyz = utils.ndgrid([x, np.r_[0], z]) + + +# plot amplitude +def plot_amplitude(ax, v): + v = spatial.vector_magnitude(v) + plt.colorbar( + ax.pcolormesh( + x, z, v.reshape(len(x), len(z), order='F'), norm=LogNorm() + ), ax=ax + ) + ax.axis('square') + ax.set_xlabel('x (m)') + ax.set_ylabel('z (m)') + + +# plot streamlines +def plot_streamlines(ax, v): + vx = v[:, 0].reshape(len(x), len(z), order='F') + vz = v[:, 2].reshape(len(x), len(z), order='F') + ax.streamplot(x, z, vx.T, vz.T, color='k') + + +# create fig, ax for electric fields and magnetic flux +fig_e, ax_e = plt.subplots( + 2, len(frequencies), figsize=(5*len(frequencies), 7) +) +fig_b, ax_b = plt.subplots( + 2, len(frequencies), figsize=(5*len(frequencies), 7) +) + +# loop over frequencies and plot +for i, frequency in enumerate(frequencies): + + # set the frequency of the dipole + edipole.frequency = frequency + + # evaluate the electric field and magnetic flux density + electric_field = edipole.electric_field(xyz) + magnetic_flux_density = edipole.magnetic_flux_density(xyz) + + # plot amplitude of electric field + for ax, reim in zip(ax_e[:, i], ['real', 'imag']): + # grab real or imag component + e_plot = getattr(electric_field, reim) + + # plot both amplitude and streamlines + plot_amplitude(ax, e_plot) + plot_streamlines(ax, e_plot) + + # set the title + ax.set_title( + 'E {} at {:1.1e} Hz'.format(reim, frequency) + ) + + # plot the amplitude of the magnetic field (note the magnetic field is into + # and out of the page in this geometry, so we don't plot vectors) + for ax, reim in zip(ax_b[:, i], ['real', 'imag']): + # grab real or imag component + b_plot = getattr(magnetic_flux_density, reim) + + # plot amplitude + plot_amplitude(ax, b_plot) + + # set the title + ax.set_title( + 'B {} at {:1.1e} Hz'.format(reim, frequency) + ) + +# format so text doesn't overlap +fig_e.tight_layout() +fig_b.tight_layout() +plt.show() + + From 32826dba6d851bf9742b53cd61b88e1e412fea57 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 3 Dec 2017 19:39:01 -0800 Subject: [PATCH 37/41] update naming in makefile to refer to geoana where it previously referred to discretize --- Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index ba895dc5..6767e42c 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ install: python setup.py install coverage: - nosetests --logging-level=INFO --with-coverage --cover-package=discretize --cover-html + nosetests --logging-level=INFO --with-coverage --cover-package=geoana --cover-html open cover/index.html lint: @@ -16,7 +16,7 @@ lint-html: pylint --output-format=html $(PACKAGE_NAME) > pylint.html graphs: - pyreverse -my -A -o pdf -p discretize discretize/**.py discretize/**/**.py + pyreverse -my -A -o pdf -p geoana geoana/**.py geoana/**/**.py tests: nosetests --logging-level=INFO From 9090eb00fd6146512812826a87faf3bddd6f5b46 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 3 Dec 2017 19:48:12 -0800 Subject: [PATCH 38/41] remove commented-out code from the static dipole --- geoana/em/static.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/geoana/em/static.py b/geoana/em/static.py index 125f26d1..346daa16 100644 --- a/geoana/em/static.py +++ b/geoana/em/static.py @@ -11,10 +11,6 @@ class MagneticDipole_WholeSpace(BaseMagneticDipole, BaseEM): - # @tr.observe('sigma') - # def _sigma_changed(self, change): - # warnings.warn("Sigma is not involved in the calculation", UserWarning) - def vector_potential(self, xyz, **kwargs): """Vector potential of a static magnetic dipole @@ -77,5 +73,3 @@ def magnetic_flux_equation(): def magnetic_field(self, xyz, **kwargs): return self.magnetic_flux(xyz, **kwargs) / self.mu - -# class From 375b8ae032b1c35a14da218b544a2bdc96b7754b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 3 Dec 2017 19:50:32 -0800 Subject: [PATCH 39/41] remove requests from setup.py --- setup.py | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.py b/setup.py index 5c30b42c..9e313113 100644 --- a/setup.py +++ b/setup.py @@ -37,7 +37,6 @@ 'numpy>=1.7', 'scipy>=0.13', 'matplotlib', - 'requests', 'ipywidgets', 'properties[math]', 'jupyter', From 9755a7e26504b8320740e3693fef1ed89eb4ca5a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 4 Dec 2017 09:59:30 -0800 Subject: [PATCH 40/41] cleanup commented out code in tests, update README --- README.rst | 6 +-- tests/test_docs.py | 16 +++---- tests/test_em_fdem_magnetic_dipole.py | 65 ++++++++++----------------- 3 files changed, 32 insertions(+), 55 deletions(-) diff --git a/README.rst b/README.rst index 16a10983..b480a994 100644 --- a/README.rst +++ b/README.rst @@ -15,14 +15,10 @@ geoana :target: https://travis-ci.org/simpeg/geoana :alt: Travis status -.. image:: https://www.quantifiedcode.com/api/v1/project/cb381a23f09245bd855f86eac295d8ec/badge.svg - :target: https://www.quantifiedcode.com/app/project/cb381a23f09245bd855f86eac295d8ec - :alt: Code issues - .. image:: https://api.codacy.com/project/badge/Grade/2e32cd28f4424dc1800f1590a64c244f :target: https://www.codacy.com/app/lindseyheagy/geoana?utm_source=github.com&utm_medium=referral&utm_content=simpeg/geoana&utm_campaign=Badge_Grade :alt: codacy status -Interactive geoscience (mostly) analytic functions. +Interactive geoscience (mostly) analytic functions in geophysics. diff --git a/tests/test_docs.py b/tests/test_docs.py index 01ab9b1e..fe31514c 100644 --- a/tests/test_docs.py +++ b/tests/test_docs.py @@ -43,14 +43,14 @@ def test_html(self): ) assert check == 0 - # def test_linkcheck(self): - # check = subprocess.call( - # ['sphinx-build', '-nW', '-b', 'linkcheck', '-d', - # '{}'.format(self.doctrees_dir), - # '{}'.format(self.docs_dir), - # '{}'.format(self.build_dir)] - # ) - # assert check == 0 + def test_linkcheck(self): + check = subprocess.call( + ['sphinx-build', '-nW', '-b', 'linkcheck', '-d', + '{}'.format(self.doctrees_dir), + '{}'.format(self.docs_dir), + '{}'.format(self.build_dir)] + ) + assert check == 0 if __name__ == '__main__': unittest.main() diff --git a/tests/test_em_fdem_magnetic_dipole.py b/tests/test_em_fdem_magnetic_dipole.py index 79e65e89..c8b0de7e 100644 --- a/tests/test_em_fdem_magnetic_dipole.py +++ b/tests/test_em_fdem_magnetic_dipole.py @@ -14,17 +14,16 @@ from discretize.utils import ndgrid, asArray_N_x_Dim -def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, epsr=1., t=0.): +def H_from_MagneticDipoleWholeSpace( + XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, + epsr=1., t=0. +): - """ - Computing magnetic fields from Magnetic Dipole in a Wholespace - TODO: - Add description of parameters - """ assert current == 1 assert loopArea == 1 assert np.all(srcLoc == np.r_[0., 0., 0.]) assert kappa == 0 + mu = mu_0 * (1+kappa) epsilon = epsilon_0 * epsr m = current * loopArea @@ -69,29 +68,28 @@ def H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1. return Hx, Hy, Hz -def B_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, epsr=1., t=0.): +def B_from_MagneticDipoleWholeSpace( + XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0, + epsr=1., t=0. +): - """ - Computing magnetic flux densites from Magnetic Dipole in a Wholespace - TODO: - Add description of parameters - """ mu = mu_0 * (1+kappa) - Hx, Hy, Hz = H_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, loopArea=loopArea, orientation=orientation, kappa=kappa, epsr=epsr) + Hx, Hy, Hz = H_from_MagneticDipoleWholeSpace( + XYZ, srcLoc, sig, f, current=current, loopArea=loopArea, + orientation=orientation, kappa=kappa, epsr=epsr + ) Bx = mu * Hx By = mu * Hy Bz = mu * Hz return Bx, By, Bz -def E_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0., epsr=1., t=0.): +def E_from_MagneticDipoleWholeSpace( + XYZ, srcLoc, sig, f, current=1., loopArea=1., orientation='X', kappa=0., + epsr=1., t=0. +): - """ - Computing analytic electric fields from Magnetic Dipole in a Wholespace - TODO: - Add description of parameters - """ mu = mu_0 * (1+kappa) epsilon = epsilon_0 * epsr m = current * loopArea @@ -108,25 +106,27 @@ def E_from_MagneticDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., loopArea=1. # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) - front = ((1j * omega(f) * mu * m) / (4.* np.pi * r**2)) * (1j * k * r + 1) * np.exp(-1j*k*r) + front = ( + ((1j * omega(f) * mu * m) / (4.* np.pi * r**2)) * + (1j * k * r + 1) * np.exp(-1j*k*r) + ) if orientation.upper() == 'X': Ey = front * (dz / r) Ez = front * (-dy / r) Ex = np.zeros_like(Ey) - return Ex, Ey, Ez elif orientation.upper() == 'Y': Ex = front * (-dz / r) Ez = front * (dx / r) Ey = np.zeros_like(Ex) - return Ex, Ey, Ez elif orientation.upper() == 'Z': Ex = front * (dy / r) Ey = front * (-dx / r) Ez = np.zeros_like(Ex) - return Ex, Ey, Ez + + return Ex, Ey, Ez class TestFDEMdipole(unittest.TestCase): @@ -349,25 +349,6 @@ class TestFDEMdipole_SimPEG(unittest.TestCase): tol = 1e-1 # error must be an order of magnitude smaller than results - # put the source at the center - - # def getFaceSrc(self, mesh): - # csx = mesh.hx.min() - # csz = mesh.hz.min() - # srcInd = ( - # (mesh.gridFz[:, 0] < csx) & - # (mesh.gridFz[:, 2] < csz/2.) & (mesh.gridFz[:, 2] > -csz/2.) - # ) - - # assert srcInd. - - # src_vecz = np.zeros(mesh.nFz, dtype=complex) - # src_vecz[srcInd] = 1. - - # return np.hstack( - # [np.zeros(mesh.vnF[:2].sum(), dtype=complex), src_vecz] - # ) - def getProjections(self, mesh): ignore_inside_radius = 10*mesh.hx.min() ignore_outside_radius = 40*mesh.hx.min() From 613767c9f747e6fe098306d672256455e51d2d4e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 12 Dec 2017 11:05:25 -0800 Subject: [PATCH 41/41] code style cleanup --- geoana/em/fdem.py | 6 +++--- geoana/em/tdem.py | 4 ++-- geoana/utils.py | 4 ++-- setup.py | 3 ++- 4 files changed, 9 insertions(+), 8 deletions(-) diff --git a/geoana/em/fdem.py b/geoana/em/fdem.py index ae958cf5..fcd04fa2 100644 --- a/geoana/em/fdem.py +++ b/geoana/em/fdem.py @@ -3,7 +3,7 @@ from __future__ import print_function from __future__ import unicode_literals -from scipy.constants import mu_0, pi, epsilon_0 +from scipy.constants import mu_0, epsilon_0 import numpy as np import warnings import properties @@ -72,8 +72,8 @@ def skin_depth(frequency, sigma, mu=mu_0): :param float sigma: electrical conductivity (S/m) :param float mu: magnetic permeability (H/m). Default: :math:`\mu_0 = 4\pi \times 10^{-7}` H/m """ - omega = omega(frequency) - return np.sqrt(2./(omega*sigma*mu)) + w = omega(frequency) + return np.sqrt(2./(w*sigma*mu)) def sigma_hat(frequency, sigma, epsilon=epsilon_0, quasistatic=False): diff --git a/geoana/em/tdem.py b/geoana/em/tdem.py index a5c78460..fa4a2a3a 100644 --- a/geoana/em/tdem.py +++ b/geoana/em/tdem.py @@ -3,7 +3,7 @@ from __future__ import print_function from __future__ import unicode_literals -from scipy.constants import mu_0, pi +from scipy.constants import mu_0 from scipy.special import erf import numpy as np import properties @@ -141,7 +141,7 @@ def magnetic_field(self, xyz): Magnetic field from an electric dipole """ dxyz = self.vector_distance(xyz) - r = self.distance(xyz) + r = self.distance(dxyz) r = spatial.repeat_scalar(r) thetar = self.theta * r diff --git a/geoana/utils.py b/geoana/utils.py index f9e69c95..83d0abf8 100644 --- a/geoana/utils.py +++ b/geoana/utils.py @@ -18,7 +18,7 @@ def mkvc(x, numDims=1): > (3, 1, 1) """ - if type(x) == np.matrix: + if isinstance(x, np.matrix): x = np.array(x) if hasattr(x, 'tovec'): @@ -71,7 +71,7 @@ def ndgrid(*args, **kwargs): # Read the keyword arguments, and only accept a vector=True/False vector = kwargs.pop('vector', True) - assert type(vector) == bool, "'vector' keyword must be a bool" + assert isinstance(vector, bool), "'vector' keyword must be a bool" assert len(kwargs) == 0, "Only 'vector' keyword accepted" # you can either pass a list [x1, x2, x3] or each seperately diff --git a/setup.py b/setup.py index 9e313113..81c3493c 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,7 @@ #!/usr/bin/env python from __future__ import print_function -"""geoana +""" +geoana Interactive geoscience (mostly) analytic functions. """