From 0838d21bce05aa79636f62441a84b7375db58604 Mon Sep 17 00:00:00 2001
From: seogi <sgkang09@gmail.com>
Date: Thu, 5 Jun 2014 11:53:15 -0700
Subject: [PATCH 1/3] Working CircularLoop code

: Can be used to compute initial condition for TDEM or
  Secondary field approach for FDEM
---
 simpegEM/Utils/Sources/CircularLoop.py | 90 ++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 simpegEM/Utils/Sources/CircularLoop.py

diff --git a/simpegEM/Utils/Sources/CircularLoop.py b/simpegEM/Utils/Sources/CircularLoop.py
new file mode 100644
index 00000000..338bc732
--- /dev/null
+++ b/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()
+
+

From 3e35c6068b315de939de362ae5fe552ddf5e3c8a Mon Sep 17 00:00:00 2001
From: seogi <sgkang09@gmail.com>
Date: Thu, 5 Jun 2014 13:11:24 -0700
Subject: [PATCH 2/3] Put CircularLoop as an option in TDEM

---
 simpegEM/TDEM/SurveyTDEM.py        | 21 +++++++++++++++++++--
 simpegEM/Utils/Sources/__init__.py |  1 -
 simpegEM/Utils/__init__.py         |  1 +
 simpegEM/__init__.py               |  2 +-
 4 files changed, 21 insertions(+), 4 deletions(-)
 delete mode 100644 simpegEM/Utils/Sources/__init__.py

diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py
index b4cc1b5d..629658de 100644
--- a/simpegEM/TDEM/SurveyTDEM.py
+++ b/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):
diff --git a/simpegEM/Utils/Sources/__init__.py b/simpegEM/Utils/Sources/__init__.py
deleted file mode 100644
index 9f44582f..00000000
--- a/simpegEM/Utils/Sources/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-from magneticDipole import MagneticDipoleVectorPotential
\ No newline at end of file
diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py
index 02b7819c..21404ecc 100644
--- a/simpegEM/Utils/__init__.py
+++ b/simpegEM/Utils/__init__.py
@@ -1,2 +1,3 @@
 import Sources
 import Ana
+import Solver
\ No newline at end of file
diff --git a/simpegEM/__init__.py b/simpegEM/__init__.py
index 5beb775b..92c6ae92 100644
--- a/simpegEM/__init__.py
+++ b/simpegEM/__init__.py
@@ -2,4 +2,4 @@
 import Utils
 import TDEM
 import FDEM
-import Base
+import Base
\ No newline at end of file

From 1bcf58154f8327a0c93a0b7059b2bcbd3fcf6c73 Mon Sep 17 00:00:00 2001
From: seogi <sgkang09@gmail.com>
Date: Thu, 5 Jun 2014 13:36:37 -0700
Subject: [PATCH 3/3] Fix import stuff

---
 simpegEM/Utils/Sources/__init__.py | 2 ++
 simpegEM/Utils/__init__.py         | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)
 create mode 100644 simpegEM/Utils/Sources/__init__.py

diff --git a/simpegEM/Utils/Sources/__init__.py b/simpegEM/Utils/Sources/__init__.py
new file mode 100644
index 00000000..6e2fedc6
--- /dev/null
+++ b/simpegEM/Utils/Sources/__init__.py
@@ -0,0 +1,2 @@
+from magneticDipole import MagneticDipoleVectorPotential
+from CircularLoop import MagneticLoopVectorPotential
diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py
index 21404ecc..c9e5e23b 100644
--- a/simpegEM/Utils/__init__.py
+++ b/simpegEM/Utils/__init__.py
@@ -1,3 +1,3 @@
 import Sources
 import Ana
-import Solver
\ No newline at end of file
+# import Solver
\ No newline at end of file