Merge branch 'master' of https://github.com/simpeg/simpegem into innerProductUpdates

simpegEM/TDEM/SurveyTDEM.py

@@ -67,7 +67,8 @@ class RxTDEM(Survey.BaseTimeRx):
 
 class TxTDEM(Survey.BaseTx):
     rxPair = RxTDEM
-    knownTxTypes = ['VMD_MVP']
+    radius = None
+    knownTxTypes = ['VMD_MVP', 'CircularLoop_MVP']
 
     def getInitialFields(self, mesh):
         F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh)
@@ -90,6 +91,23 @@ class TxTDEM(Survey.BaseTx):
 
         return {"b": mesh.edgeCurl*MVP}
 
+    def _getInitialFields_CircularLoop_MVP(self, mesh):
+        """Circular Loop, magnetic vector potential"""
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = Sources.MagneticLoopVectorPotential(self.loc, mesh.gridEy, 'y', self.radius)
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVPx = Sources.MagneticLoopVectorPotential(self.loc, mesh.gridEx, 'x', self.radius)
+            MVPy = Sources.MagneticLoopVectorPotential(self.loc, mesh.gridEy, 'y', self.radius)
+            MVPz = Sources.MagneticLoopVectorPotential(self.loc, mesh.gridEz, 'z', self.radius)
+            MVP = np.concatenate((MVPx, MVPy, MVPz))
+        else:
+            raise Exception('Unknown mesh for CircularLoop')
+
+        return {"b": mesh.edgeCurl*MVP}        
+
     def getJs(self, mesh, time):
         return None
 
@@ -97,7 +115,6 @@ class SurveyTDEM(Survey.BaseSurvey):
     """
         docstring for SurveyTDEM
     """
-
     txPair = TxTDEM
 
     def __init__(self, txList, **kwargs):

simpegEM/Utils/Sources/CircularLoop.py

@@ -0,0 +1,90 @@
+from SimPEG import *
+from scipy.special import ellipk, ellipe
+
+def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius):
+    """
+        Calculate the vector potential of horizontal circular loop
+        at given locations
+
+        :param numpy.ndarray txLoc: Location of the transmitter(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 I: Input current of the loop
+        :param numpy.ndarray radius: radius of the loop
+        :rtype: numpy.ndarray
+        :return: The vector potential each dipole at each observation location
+    """
+
+    txLoc = np.atleast_2d(txLoc)
+    obsLoc = np.atleast_2d(obsLoc)
+
+    n = obsLoc.shape[0]
+    nTx = txLoc.shape[0]
+
+    if component=='z':
+        A = np.zeros((n, nTx))
+        if nTx ==1:
+            return A.flatten()        
+        return A
+
+    else:
+
+        A = np.zeros((n, nTx))
+        for i in range (nTx):
+            x = obsLoc[:, 0] - txLoc[i, 0]
+            y = obsLoc[:, 1] - txLoc[i, 1]
+            z = obsLoc[:, 2] - txLoc[i, 2]
+            r = np.sqrt(x**2 + y**2)
+            m = (4 * radius * r) / ((radius + r)**2 + z**2)
+            m[m > 1.] = 1.
+            # m might be slightly larger than 1 due to rounding errors
+            # but ellipke requires 0 <= m <= 1
+            K = ellipk(m)
+            E = ellipe(m)
+            ind = (r > 0) & (m < 1)
+            # % 1/r singular at r = 0 and K(m) singular at m = 1
+            Aphi = np.zeros(n)
+            # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
+            Aphi[ind] = 4e-7 / np.sqrt(m[ind])  * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
+            if component == 'x':
+                A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
+            elif component == 'y':
+                A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
+            else:
+                raise ValueError('Invalid component')               
+
+        if nTx == 1:
+            return A.flatten()                
+        return A
+
+if __name__ == '__main__':
+    from SimPEG import Mesh
+    import matplotlib.pyplot as plt
+    cs = 20
+    ncx, ncy, ncz = 41, 41, 40
+    hx = np.ones(ncx)*cs
+    hy = np.ones(ncy)*cs
+    hz = np.ones(ncz)*cs
+    mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
+    txLoc = np.r_[0., 0., 0.]
+    Ax = MagneticLoopVectorPotential(txLoc, mesh.gridEx, 'x', 200)
+    Ay = MagneticLoopVectorPotential(txLoc, mesh.gridEy, 'y', 200)
+    Az = MagneticLoopVectorPotential(txLoc, mesh.gridEz, 'z', 200)
+    A = np.r_[Ax, Ay, Az]
+    B0 = mesh.edgeCurl*A
+    J0 = mesh.edgeCurl.T*B0
+
+    # mesh.plotImage(A, vType = 'Ex')
+    # mesh.plotImage(A, vType = 'Ey')
+
+    mesh.plotImage(B0, vType = 'Fx')
+    mesh.plotImage(B0, vType = 'Fy')
+    mesh.plotImage(B0, vType = 'Fz')
+
+    # # mesh.plotImage(J0, vType = 'Ex')
+    # mesh.plotImage(J0, vType = 'Ey')
+    # mesh.plotImage(J0, vType = 'Ez')
+
+    plt.show()
+
+

simpegEM/Utils/Sources/__init__.py

@@ -1 +1,2 @@
-from magneticDipole import MagneticDipoleVectorPotential
+from magneticDipole import MagneticDipoleVectorPotential
+from CircularLoop import MagneticLoopVectorPotential

simpegEM/Utils/__init__.py

@@ -1,2 +1,3 @@
 import Sources
 import Ana
+# import Solver

simpegEM/__init__.py

@@ -2,4 +2,4 @@
 import Utils
 import TDEM
 import FDEM
-import Base
+import Base