Merge pull request #13 from simpeg/magnetic_pole

magnetic pole

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.0.4
+current_version = 0.0.5
 files = setup.py geoana/__init__.py docs/conf.py
 

appveyor.yml

@@ -38,8 +38,8 @@ install:
   - conda info -a
   - "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
+  - pip install -r requirements_dev.txt
 
 test_script:
   - nosetests tests -v -s

docs/conf.py

@@ -64,9 +64,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'0.0.4'
+version = u'0.0.5'
 # The full version, including alpha/beta/rc tags.
-release = u'0.0.4'
+release = u'0.0.5'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

examples/em/plot_static_magnetic_fields_from_pole_source.py

@@ -0,0 +1,120 @@
+"""
+Total magnetic fields: Dipole and Pole sources
+==============================================
+
+In this example, we plot anomalous total magnetic field
+from a magnetic dipole and pole targets. These targets are
+excited by Earth magnetic fields.
+We can vary the direction of the Earth magnetic field, and
+magnetic moment of the target.
+
+:author: Seogi Kang (`@sgkang <https://github.com/sgkang>`_)
+:date: Aug 19, 2018
+
+"""
+
+
+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 static
+
+###############################################################################
+# Setup
+# -----
+#
+# define the location, orientation, and source, physical properties of the
+# wholespace and source parameters
+
+mu = mu_0  # permeability of free space (this is the default)
+location = np.r_[0., 0., -10.]  # location of the dipole or pole
+
+# dipole parameters
+moment = 1
+# inclination and declination (e.g. Vancouver)
+inclination, declination = 67., 0.
+
+
+###############################################################################
+# Magnetostatic Dipole and Loop
+# -----------------------------
+#
+# Here, we build the geoana magnetic dipole and poie in a wholespace
+# using the parameters defined above.
+
+def id_to_cartesian(inclination, declination):
+    ux = np.cos(inclination/180.*np.pi)*np.sin(declination/180.*np.pi)
+    uy = np.cos(inclination/180.*np.pi)*np.cos(declination/180.*np.pi)
+    uz = -np.sin(inclination/180.*np.pi)
+    return np.r_[ux, uy, uz]
+
+orientation = id_to_cartesian(inclination, declination)
+
+dipole = static.MagneticDipoleWholeSpace(
+    location=location,
+    orientation=orientation,
+    moment=moment
+)
+
+pole = static.MagneticPoleWholeSpace(
+    location=location,
+    orientation=orientation,
+    moment=moment
+)
+
+###############################################################################
+# Evaluate magnetic fields
+# --------------------------
+#
+# Next, we construct a grid where we want to plot the magentic fields and
+# evaluate
+
+x = np.linspace(-36, 36, 100)
+y = np.linspace(-36, 36, 100)
+xyz = utils.ndgrid([x, y, np.r_[1.]])
+
+# evaluate the magnetic field
+b_vec_dipole = dipole.magnetic_flux_density(xyz)
+b_vec_pole = pole.magnetic_flux_density(xyz)
+b_total_dipole = dipole.dot_orientation(b_vec_dipole)
+b_total_pole = pole.dot_orientation(b_vec_pole)
+###############################################################################
+#
+# and define plotting code to plot an image of the amplitude of the vector
+# field / flux as well as the streamlines
+
+
+def plot_amplitude(ax, v):
+    plt.colorbar(
+        ax.pcolormesh(
+            x, y, v.reshape(len(x), len(y), order='F')
+        ), ax=ax
+    )
+    ax.axis('square')
+    ax.set_xlabel('y (east,  m)')
+    ax.set_ylabel('x (north,  m)')
+
+###############################################################################
+#
+# Create subplots for plotting the results. Loop over frequencies and plot the
+# electric and magnetic fields along a slice through the center of the dipole.
+
+fig, ax = plt.subplots(1, 2, figsize=(12, 5))
+
+# plot dipole vector potential
+plot_amplitude(ax[0], b_total_dipole)
+
+# plot loop vector potential
+plot_amplitude(ax[1], b_total_pole)
+
+
+# set the titles
+ax[0].set_title("Total field: dipole")
+ax[1].set_title("Total field: pole")
+
+# format so text doesn't overlap
+plt.tight_layout()
+plt.show()

geoana/__init__.py

@@ -2,7 +2,7 @@
 from . import em
 from . import utils
 
-__version__ = '0.0.4'
+__version__ = '0.0.5'
 __author__ = 'SimPEG developers'
 __license__ = 'MIT'
 __copyright__ = 'Copyright 2017-2018 SimPEG developers'

geoana/em/base.py

@@ -83,8 +83,8 @@ def cross_orientation(self, xyz):
         """
         orientation = np.kron(
             np.atleast_2d(
-                np.array(self.orientation)), np.ones((xyz.shape[0], 1)
-            )
+                np.array(self.orientation)
+            ), np.ones((xyz.shape[0], 1))
         )
         return np.cross(xyz, orientation)
 

geoana/em/static.py

@@ -12,7 +12,10 @@
 from .base import BaseEM, BaseDipole, BaseMagneticDipole, BaseElectricDipole
 from .. import spatial
 
-__all__ = ["MagneticDipoleWholeSpace", "CircularLoopWholeSpace"]
+__all__ = [
+    "MagneticDipoleWholeSpace", "CircularLoopWholeSpace",
+    "MagneticPoleWholeSpace"
+]
 
 
 class MagneticDipoleWholeSpace(BaseMagneticDipole, BaseEM):
@@ -73,7 +76,6 @@ def vector_potential(self, xyz, coordinates="cartesian"):
 
         return a
 
-
     def magnetic_flux_density(self, xyz, coordinates="cartesian"):
         """Magnetic flux (:math:`\\vec{b}`) of a static magnetic dipole
 
@@ -152,6 +154,77 @@ def magnetic_field(self, xyz, coordinates="cartesian"):
         return self.magnetic_flux(xyz, coordinates=coordinates) / self.mu
 
 
+class MagneticPoleWholeSpace(BaseMagneticDipole, BaseEM):
+    """
+    Static magnetic pole in a wholespace.
+    """
+
+    def magnetic_flux_density(self, xyz, coordinates="cartesian"):
+        """Magnetic flux (:math:`\\vec{b}`) of a static magnetic dipole
+
+        **Required**
+
+        :param numpy.ndarray xyz: Location of the receivers(s)
+
+        **Optional**
+
+        :param str coordinates: coordinate system that the xyz is provided
+                                in and that the solution will be returned
+                                in (cartesian or cylindrical).
+                                Default: `"cartesian"`
+
+        **Returns**
+
+        :rtype: numpy.ndarray
+        :return: The magnetic flux at each observation location
+        """
+
+        supported_coordinates = ["cartesian", "cylindrical"]
+        assert coordinates.lower() in supported_coordinates, (
+            "coordinates must be in {}, the coordinate system "
+            "you provided, {}, is not yet supported".format(
+                supported_coordinates, coordinates
+            )
+        )
+
+        n_obs = xyz.shape[0]
+
+        if coordinates.lower() == "cylindrical":
+            xyz = spatial.cylindrical_2_cartesian(xyz)
+
+        r = self.vector_distance(xyz)
+        dxyz = spatial.repeat_scalar(self.distance(xyz))
+
+        b = self.moment * self.mu / (4 * np.pi * (dxyz ** 3)) * r
+
+        if coordinates.lower() == "cylindrical":
+            b = spatial.cartesian_2_cylindrical(xyz, b)
+
+        return b
+
+    def magnetic_field(self, xyz, coordinates="cartesian"):
+        """Magnetic field (:math:`\\vec{h}`) of a static magnetic dipole
+
+        **Required**
+
+        :param numpy.ndarray xyz: Location of the receivers(s)
+
+        **Optional**
+
+        :param str coordinates: coordinate system that the xyz is provided
+                                in and that the solution will be returned
+                                in (cartesian or cylindrical).
+                                Default: `"cartesian"`
+
+        **Returns**
+
+        :rtype: numpy.ndarray
+        :return: The magnetic field at each observation location
+
+        """
+        return self.magnetic_flux(xyz, coordinates=coordinates) / self.mu
+
+
 class CircularLoopWholeSpace(BaseDipole, BaseEM):
 
     """

geoana/spatial.py

@@ -52,6 +52,7 @@ def cylindrical_2_cartesian(grid, vec=None):
 
     return np.vstack(newvec).T
 
+
 def cartesian_2_cylindrical(grid, vec=None):
     """
     Takes a grid or vector (if provided)

setup.py

@@ -31,7 +31,7 @@
 
 setup(
     name = 'geoana',
-    version = '0.0.4',
+    version = '0.0.5',
     packages = find_packages(),
     install_requires = [
         'future',