Merge branch 'master' of https://github.com/simpeg/simpeg into mesh/tree

.travis.yml

@@ -5,7 +5,11 @@ python:
 sudo: false
 
 env:
-    - TEST_DIR=tests/em
+    - TEST_DIR=tests/em/examples
+    - TEST_DIR=tests/em/fdem/forward
+    - TEST_DIR=tests/em/fdem/inverse/derivs
+    - TEST_DIR=tests/em/fdem/inverse/adjoint
+    - TEST_DIR=tests/em/tdem
     - TEST_DIR=tests/mesh
     - TEST_DIR=tests/flow
     - TEST_DIR=tests/utils

SimPEG/EM/FDEM/FDEM.py

@@ -1,9 +1,9 @@
-from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
+from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
-from SurveyFDEM import SurveyFDEM
-from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j
+from SurveyFDEM import Survey as SurveyFDEM
+from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j
 from SimPEG.EM.Base import BaseEMProblem
-from SimPEG.EM.Utils.EMUtils import omega
+from SimPEG.EM.Utils import omega
 
 
 class BaseFDEMProblem(BaseEMProblem):
@@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
             {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`ProblemFDEM_e`
-        or :code:`ProblemFDEM_b`) or the magnetic field
+        if using the E-B formulation (:code:`Problem_e`
+        or :code:`Problem_b`) or the magnetic field
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
         .. math ::
@@ -26,13 +26,13 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`).
+        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`).
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
 
     surveyPair = SurveyFDEM
-    fieldsPair = FieldsFDEM
+    fieldsPair = Fields
 
     def fields(self, m=None):
         """
@@ -185,7 +185,7 @@ class BaseFDEMProblem(BaseEMProblem):
 ################################ E-B Formulation #########################################
 ##########################################################################################
 
-class ProblemFDEM_e(BaseFDEMProblem):
+class Problem_e(BaseFDEMProblem):
     """
         By eliminating the magnetic flux density using
 
@@ -205,7 +205,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
 
     _fieldType = 'e'
     _eqLocs    = 'FE'
-    fieldsPair = FieldsFDEM_e
+    fieldsPair = Fields_e
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -284,7 +284,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
                 return None
 
 
-class ProblemFDEM_b(BaseFDEMProblem):
+class Problem_b(BaseFDEMProblem):
     """
         We eliminate \\\(\\\mathbf{e}\\\) using
 
@@ -304,7 +304,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
 
     _fieldType = 'b'
     _eqLocs    = 'FE'
-    fieldsPair = FieldsFDEM_b
+    fieldsPair = Fields_b
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -425,7 +425,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
 ##########################################################################################
 
 
-class ProblemFDEM_j(BaseFDEMProblem):
+class Problem_j(BaseFDEMProblem):
     """
         We eliminate \\\(\\\mathbf{h}\\\) using
 
@@ -446,7 +446,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
 
     _fieldType = 'j'
     _eqLocs    = 'EF'
-    fieldsPair = FieldsFDEM_j
+    fieldsPair = Fields_j
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -558,7 +558,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
 
 
 
-class ProblemFDEM_h(BaseFDEMProblem):
+class Problem_h(BaseFDEMProblem):
     """
         We eliminate \\\(\\\mathbf{j}\\\) using
 
@@ -576,7 +576,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
 
     _fieldType = 'h'
     _eqLocs    = 'EF'
-    fieldsPair = FieldsFDEM_h
+    fieldsPair = Fields_h
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -1,13 +1,16 @@
-from SimPEG import Survey, Problem, Utils, np, sp
-from SimPEG.EM.Utils.EMUtils import omega
+import numpy as np
+import scipy.sparse as sp
+import SimPEG
+from SimPEG import Utils
+from SimPEG.EM.Utils import omega
 
 
-class FieldsFDEM(Problem.Fields):
+class Fields(SimPEG.Problem.Fields):
     """Fancy Field Storage for a FDEM survey."""
     knownFields = {}
     dtype = complex
 
-class FieldsFDEM_e(FieldsFDEM):
+class Fields_e(Fields):
     knownFields = {'eSolution':'E'}
     aliasFields = {
                     'e' : ['eSolution','E','_e'],
@@ -19,7 +22,7 @@ class FieldsFDEM_e(FieldsFDEM):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
@@ -89,7 +92,7 @@ class FieldsFDEM_e(FieldsFDEM):
         return self._bSecondaryDeriv_m(src, v, adjoint)
 
 
-class FieldsFDEM_b(FieldsFDEM):
+class Fields_b(Fields):
     knownFields = {'bSolution':'F'}
     aliasFields = {
                     'b' : ['bSolution','F','_b'],
@@ -101,7 +104,7 @@ class FieldsFDEM_b(FieldsFDEM):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
@@ -190,7 +193,7 @@ class FieldsFDEM_b(FieldsFDEM):
         return self._eSecondaryDeriv_m(src, v, adjoint)
 
 
-class FieldsFDEM_j(FieldsFDEM):
+class Fields_j(Fields):
     knownFields = {'jSolution':'F'}
     aliasFields = {
                     'j' : ['jSolution','F','_j'],
@@ -202,7 +205,7 @@ class FieldsFDEM_j(FieldsFDEM):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
@@ -293,7 +296,7 @@ class FieldsFDEM_j(FieldsFDEM):
         return self._hSecondaryDeriv_m(src, v, adjoint)
 
 
-class FieldsFDEM_h(FieldsFDEM):
+class Fields_h(Fields):
     knownFields = {'hSolution':'E'}
     aliasFields = {
                     'h' : ['hSolution','E','_h'],
@@ -305,7 +308,7 @@ class FieldsFDEM_h(FieldsFDEM):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob

SimPEG/EM/FDEM/SrcFDEM.py

@@ -0,0 +1,347 @@
+from SimPEG import Survey, Problem, Utils, np, sp
+from scipy.constants import mu_0
+from SimPEG.EM.Utils import *
+# from SurveyFDEM import Rx
+ 
+
+class BaseSrc(Survey.BaseSrc):
+    freq = None
+    # rxPair = Rx
+    integrate = True
+
+    def eval(self, prob):
+        S_m = self.S_m(prob)
+        S_e = self.S_e(prob)
+        return S_m, S_e
+
+    def evalDeriv(self, prob, v, adjoint=False):
+        return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
+
+    def bPrimary(self, prob):
+        return None
+
+    def hPrimary(self, prob):
+        return None
+
+    def ePrimary(self, prob):
+        return None
+
+    def jPrimary(self, prob):
+        return None
+
+    def S_m(self, prob):
+        return None
+
+    def S_e(self, prob):
+        return None
+
+    def S_mDeriv(self, prob, v, adjoint = False):
+        return None
+
+    def S_eDeriv(self, prob, v, adjoint = False):
+        return None
+
+
+class RawVec_e(BaseSrc):
+    """
+        RawVec electric source. It is defined by the user provided vector S_e
+
+        :param numpy.array S_e: electric source term
+        :param float freq: frequency
+        :param rxList: receiver list
+    """
+
+    def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+        self._S_e = np.array(S_e,dtype=complex)
+        self._ePrimary = ePrimary
+        self._bPrimary = bPrimary
+        self._hPrimary = hPrimary
+        self._jPrimary = jPrimary
+        self.freq = float(freq)
+        BaseSrc.__init__(self, rxList)
+
+    def S_e(self, prob):
+        return self._S_e
+
+    def ePrimary(self, prob):
+        return self._ePrimary
+
+    def bPrimary(self, prob):
+        return self._bPrimary
+
+    def hPrimary(self, prob):
+        return self._hPrimary
+
+    def jPrimary(self, prob):
+        return self._jPrimary
+
+
+class RawVec_m(BaseSrc):
+    """
+        RawVec magnetic source. It is defined by the user provided vector S_m
+
+        :param numpy.array S_m: magnetic source term
+        :param float freq: frequency
+        :param rxList: receiver list
+    """
+
+    def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+        self._S_m = np.array(S_m,dtype=complex)
+        self.freq = float(freq)
+        self.integrate = integrate
+        self._ePrimary = np.array(ePrimary,dtype=complex)
+        self._bPrimary = np.array(bPrimary,dtype=complex)
+        self._hPrimary = np.array(hPrimary,dtype=complex)
+        self._jPrimary = np.array(jPrimary,dtype=complex)
+
+        BaseSrc.__init__(self, rxList)
+
+    def S_m(self, prob):
+        return self._S_m
+
+    def ePrimary(self, prob):
+        return self._ePrimary
+
+    def bPrimary(self, prob):
+        return self._bPrimary
+
+    def hPrimary(self, prob):
+        return self._hPrimary
+
+    def jPrimary(self, prob):
+        return self._jPrimary
+
+
+class RawVec(BaseSrc):
+    """
+        RawVec source. It is defined by the user provided vectors S_m, S_e
+
+        :param numpy.array S_m: magnetic source term
+        :param numpy.array S_e: electric source term
+        :param float freq: frequency
+        :param rxList: receiver list
+    """
+    def __init__(self, rxList, freq, S_m, S_e, integrate = True):
+        self._S_m = np.array(S_m,dtype=complex)
+        self._S_e = np.array(S_e,dtype=complex)
+        self.freq = float(freq)
+        self.integrate = integrate
+        BaseSrc.__init__(self, rxList)
+
+    def S_m(self, prob):
+        if prob._eqLocs is 'EF' and self.integrate is True:
+            return prob.Me * self._S_m
+        return self._S_m
+
+    def S_e(self, prob):
+        if prob._eqLocs is 'FE' and self.integrate is True:
+            return prob.Me * self._S_e
+        return self._S_e
+
+
+class MagDipole(BaseSrc):
+
+    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
+        self.freq = float(freq)
+        self.loc = loc
+        self.orientation = orientation
+        self.moment = moment
+        self.mu = mu
+        self.integrate = False
+        BaseSrc.__init__(self, rxList)
+
+    def bPrimary(self, prob):
+        eqLocs = prob._eqLocs
+
+        if eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                # TODO ?
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+
+        else:
+            srcfct = MagneticDipoleVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
+    def hPrimary(self, prob):
+        b = self.bPrimary(prob)
+        return h_from_b(prob,b)
+
+    def S_m(self, prob):
+        b_p = self.bPrimary(prob)
+        return -1j*omega(self.freq)*b_p
+
+    def S_e(self, prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
+
+
+class MagDipole_Bfield(BaseSrc):
+
+    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
+    #TODO: neither does moment
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
+        self.freq = float(freq)
+        self.loc = loc
+        self.orientation = orientation
+        self.moment = moment
+        self.mu = mu
+        BaseSrc.__init__(self, rxList)
+
+    def bPrimary(self, prob):
+        eqLocs = prob._eqLocs
+
+        if eqLocs is 'FE':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl
+
+        elif eqLocs is 'EF':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl.T
+
+        srcfct = MagneticDipoleFields
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                # TODO ?
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            b = np.concatenate((bx,bz))
+        else:
+            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            b = np.concatenate((bx,by,bz))
+
+        return b
+
+    def hPrimary(self, prob):
+        b = self.bPrimary(prob)
+        return h_from_b(prob, b)
+
+    def S_m(self, prob):
+        b = self.bPrimary(prob)
+        return -1j*omega(self.freq)*b
+
+    def S_e(self, prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
+
+
+class CircularLoop(BaseSrc):
+
+    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
+    def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
+        self.freq = float(freq)
+        self.orientation = orientation
+        self.radius = radius
+        self.mu = mu
+        self.loc = loc
+        self.integrate = False
+        BaseSrc.__init__(self, rxList)
+
+    def bPrimary(self, prob):
+        eqLocs = prob._eqLocs
+
+        if eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                # TODO ?
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+
+        else:
+            srcfct = MagneticDipoleVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
+            ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
+            az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
+    def hPrimary(self, prob):
+        b = self.bPrimary(prob)
+        return 1./self.mu*b
+
+    def S_m(self, prob):
+        b = self.bPrimary(prob)
+        return -1j*omega(self.freq)*b
+
+    def S_e(self, prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
+
+

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -1,13 +1,13 @@
-from SimPEG import Survey, Problem, Utils, np, sp
-from SimPEG.EM.Utils import SrcUtils
-from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
+import SimPEG
+from SimPEG.EM.Utils import *
 from scipy.constants import mu_0
+import SrcFDEM as Src
 
 ####################################################
 # Receivers
 ####################################################
 
-class RxFDEM(Survey.BaseRx):
+class Rx(SimPEG.Survey.BaseRx):
 
     knownRxTypes = {
                     'exr':['e', 'Ex', 'real'],
@@ -41,7 +41,7 @@ class RxFDEM(Survey.BaseRx):
     radius = None
 
     def __init__(self, locs, rxType):
-        Survey.BaseRx.__init__(self, locs, rxType)
+        SimPEG.Survey.BaseRx.__init__(self, locs, rxType)
 
     @property
     def projField(self):
@@ -87,366 +87,21 @@ class RxFDEM(Survey.BaseRx):
         return Pv
 
 
-####################################################
-# Sources
-####################################################
-
-class SrcFDEM(Survey.BaseSrc):
-    freq = None
-    rxPair = RxFDEM
-    integrate = True
-
-    def eval(self, prob):
-        S_m = self.S_m(prob)
-        S_e = self.S_e(prob)
-        return S_m, S_e
-
-    def evalDeriv(self, prob, v, adjoint=False):
-        return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
-
-    def bPrimary(self, prob):
-        return None
-
-    def hPrimary(self, prob):
-        return None
-
-    def ePrimary(self, prob):
-        return None
-
-    def jPrimary(self, prob):
-        return None
-
-    def S_m(self, prob):
-        return None
-
-    def S_e(self, prob):
-        return None
-
-    def S_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def S_eDeriv(self, prob, v, adjoint = False):
-        return None
-
-
-class SrcFDEM_RawVec_e(SrcFDEM):
-    """
-        RawVec electric source. It is defined by the user provided vector S_e
-
-        :param numpy.array S_e: electric source term
-        :param float freq: frequency
-        :param rxList: receiver list
-    """
-
-    def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
-        self._S_e = np.array(S_e,dtype=complex)
-        self._ePrimary = ePrimary
-        self._bPrimary = bPrimary
-        self._hPrimary = hPrimary
-        self._jPrimary = jPrimary
-        self.freq = float(freq)
-        SrcFDEM.__init__(self, rxList)
-
-    def S_e(self, prob):
-        return self._S_e
-
-    def ePrimary(self, prob):
-        return self._ePrimary
-
-    def bPrimary(self, prob):
-        return self._bPrimary
-
-    def hPrimary(self, prob):
-        return self._hPrimary
-
-    def jPrimary(self, prob):
-        return self._jPrimary
-
-
-class SrcFDEM_RawVec_m(SrcFDEM):
-    """
-        RawVec magnetic source. It is defined by the user provided vector S_m
-
-        :param numpy.array S_m: magnetic source term
-        :param float freq: frequency
-        :param rxList: receiver list
-    """
-
-    def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
-        self._S_m = np.array(S_m,dtype=complex)
-        self.freq = float(freq)
-        self.integrate = integrate
-        self._ePrimary = np.array(ePrimary,dtype=complex)
-        self._bPrimary = np.array(bPrimary,dtype=complex)
-        self._hPrimary = np.array(hPrimary,dtype=complex)
-        self._jPrimary = np.array(jPrimary,dtype=complex)
-
-        SrcFDEM.__init__(self, rxList)
-
-    def S_m(self, prob):
-        return self._S_m
-
-    def ePrimary(self, prob):
-        return self._ePrimary
-
-    def bPrimary(self, prob):
-        return self._bPrimary
-
-    def hPrimary(self, prob):
-        return self._hPrimary
-
-    def jPrimary(self, prob):
-        return self._jPrimary
-
-
-class SrcFDEM_RawVec(SrcFDEM):
-    """
-        RawVec source. It is defined by the user provided vectors S_m, S_e
-
-        :param numpy.array S_m: magnetic source term
-        :param numpy.array S_e: electric source term
-        :param float freq: frequency
-        :param rxList: receiver list
-    """
-    def __init__(self, rxList, freq, S_m, S_e, integrate = True):
-        self._S_m = np.array(S_m,dtype=complex)
-        self._S_e = np.array(S_e,dtype=complex)
-        self.freq = float(freq)
-        self.integrate = integrate
-        SrcFDEM.__init__(self, rxList)
-
-    def S_m(self, prob):
-        if prob._eqLocs is 'EF' and self.integrate is True:
-            return prob.Me * self._S_m
-        return self._S_m
-
-    def S_e(self, prob):
-        if prob._eqLocs is 'FE' and self.integrate is True:
-            return prob.Me * self._S_e
-        return self._S_e
-
-
-class SrcFDEM_MagDipole(SrcFDEM):
-
-    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
-        self.freq = float(freq)
-        self.loc = loc
-        self.orientation = orientation
-        self.moment = moment
-        self.mu = mu
-        self.integrate = False
-        SrcFDEM.__init__(self, rxList)
-
-    def bPrimary(self, prob):
-        eqLocs = prob._eqLocs
-
-        if eqLocs is 'FE':
-            gridX = prob.mesh.gridEx
-            gridY = prob.mesh.gridEy
-            gridZ = prob.mesh.gridEz
-            C = prob.mesh.edgeCurl
-
-        elif eqLocs is 'EF':
-            gridX = prob.mesh.gridFx
-            gridY = prob.mesh.gridFy
-            gridZ = prob.mesh.gridFz
-            C = prob.mesh.edgeCurl.T
-
-
-        if prob.mesh._meshType is 'CYL':
-            if not prob.mesh.isSymmetric:
-                # TODO ?
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
-
-        else:
-            srcfct = SrcUtils.MagneticDipoleVectorPotential
-            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
-            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
-            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
-            a = np.concatenate((ax, ay, az))
-
-        return C*a
-
-    def hPrimary(self, prob):
-        b = self.bPrimary(prob)
-        return h_from_b(prob,b)
-
-    def S_m(self, prob):
-        b_p = self.bPrimary(prob)
-        return -1j*omega(self.freq)*b_p
-
-    def S_e(self, prob):
-        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
-            return None
-        else:
-            eqLocs = prob._eqLocs
-
-            if eqLocs is 'FE':
-                mui_s = prob.curModel.mui - 1./self.mu
-                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
-                C = prob.mesh.edgeCurl
-            elif eqLocs is 'EF':
-                mu_s = prob.curModel.mu - self.mu
-                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
-                C = prob.mesh.edgeCurl.T
-
-            return -C.T * (MMui_s * self.bPrimary(prob))
-
-
-class SrcFDEM_MagDipole_Bfield(SrcFDEM):
-
-    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
-    #TODO: neither does moment
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
-        self.freq = float(freq)
-        self.loc = loc
-        self.orientation = orientation
-        self.moment = moment
-        self.mu = mu
-        SrcFDEM.__init__(self, rxList)
-
-    def bPrimary(self, prob):
-        eqLocs = prob._eqLocs
-
-        if eqLocs is 'FE':
-            gridX = prob.mesh.gridFx
-            gridY = prob.mesh.gridFy
-            gridZ = prob.mesh.gridFz
-            C = prob.mesh.edgeCurl
-
-        elif eqLocs is 'EF':
-            gridX = prob.mesh.gridEx
-            gridY = prob.mesh.gridEy
-            gridZ = prob.mesh.gridEz
-            C = prob.mesh.edgeCurl.T
-
-        srcfct = SrcUtils.MagneticDipoleFields
-        if prob.mesh._meshType is 'CYL':
-            if not prob.mesh.isSymmetric:
-                # TODO ?
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
-            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
-            b = np.concatenate((bx,bz))
-        else:
-            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
-            by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
-            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
-            b = np.concatenate((bx,by,bz))
-
-        return b
-
-    def hPrimary(self, prob):
-        b = self.bPrimary(prob)
-        return h_from_b(prob, b)
-
-    def S_m(self, prob):
-        b = self.bPrimary(prob)
-        return -1j*omega(self.freq)*b
-
-    def S_e(self, prob):
-        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
-            return None
-        else:
-            eqLocs = prob._eqLocs
-
-            if eqLocs is 'FE':
-                mui_s = prob.curModel.mui - 1./self.mu
-                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
-                C = prob.mesh.edgeCurl
-            elif eqLocs is 'EF':
-                mu_s = prob.curModel.mu - self.mu
-                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
-                C = prob.mesh.edgeCurl.T
-
-            return -C.T * (MMui_s * self.bPrimary(prob))
-
-
-class SrcFDEM_CircularLoop(SrcFDEM):
-
-    #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
-    def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
-        self.freq = float(freq)
-        self.orientation = orientation
-        self.radius = radius
-        self.mu = mu
-        self.loc = loc
-        self.integrate = False
-        SrcFDEM.__init__(self, rxList)
-
-    def bPrimary(self, prob):
-        eqLocs = prob._eqLocs
-
-        if eqLocs is 'FE':
-            gridX = prob.mesh.gridEx
-            gridY = prob.mesh.gridEy
-            gridZ = prob.mesh.gridEz
-            C = prob.mesh.edgeCurl
-
-        elif eqLocs is 'EF':
-            gridX = prob.mesh.gridFx
-            gridY = prob.mesh.gridFy
-            gridZ = prob.mesh.gridFz
-            C = prob.mesh.edgeCurl.T
-
-        if prob.mesh._meshType is 'CYL':
-            if not prob.mesh.isSymmetric:
-                # TODO ?
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
-
-        else:
-            srcfct = SrcUtils.MagneticDipoleVectorPotential
-            ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
-            ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
-            az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
-            a = np.concatenate((ax, ay, az))
-
-        return C*a
-
-    def hPrimary(self, prob):
-        b = self.bPrimary(prob)
-        return 1./self.mu*b
-
-    def S_m(self, prob):
-        b = self.bPrimary(prob)
-        return -1j*omega(self.freq)*b
-
-    def S_e(self, prob):
-        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
-            return None
-        else:
-            eqLocs = prob._eqLocs
-
-            if eqLocs is 'FE':
-                mui_s = prob.curModel.mui - 1./self.mu
-                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
-                C = prob.mesh.edgeCurl
-            elif eqLocs is 'EF':
-                mu_s = prob.curModel.mu - self.mu
-                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
-                C = prob.mesh.edgeCurl.T
-
-            return -C.T * (MMui_s * self.bPrimary(prob))
-
-
 ####################################################
 # Survey
 ####################################################
 
-class SurveyFDEM(Survey.BaseSurvey):
+class Survey(SimPEG.Survey.BaseSurvey):
     """
         docstring for SurveyFDEM
     """
 
-    srcPair = SrcFDEM
+    srcPair = Src.BaseSrc
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        Survey.BaseSurvey.__init__(self, **kwargs)
+        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
         _freqDict = {}
         for src in srcList:
@@ -481,7 +136,7 @@ class SurveyFDEM(Survey.BaseSurvey):
         return self._freqDict[freq]
 
     def projectFields(self, u):
-        data = Survey.Data(self)
+        data = SimPEG.Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
                 data[src, rx] = rx.projectFields(src, self.mesh, u)

SimPEG/EM/FDEM/__init__.py

@@ -1,3 +1,3 @@
-from SurveyFDEM import *
-from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
+from SurveyFDEM import Rx, Src, Survey
+from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h
 from FieldsFDEM import *

SimPEG/EM/TDEM/BaseTDEM.py

@@ -1,6 +1,6 @@
 from SimPEG import Solver, Problem
 from SimPEG.Problem import BaseTimeProblem
-from SimPEG.EM.Utils import SrcUtils
+from SimPEG.EM.Utils import *
 from scipy.constants import mu_0
 from SimPEG.Utils import sdiag, mkvc
 from SimPEG import Utils, Mesh

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -1,6 +1,6 @@
 from SimPEG import Utils, Survey, np
 from SimPEG.Survey import BaseSurvey
-from SimPEG.EM.Utils import SrcUtils
+from SimPEG.EM.Utils import *
 from BaseTDEM import FieldsTDEM
 
 
@@ -87,11 +87,11 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         """Vertical magnetic dipole, magnetic vector potential"""
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
-                MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
             else:
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
         elif mesh._meshType is 'TENSOR':
-            MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
             raise Exception('Unknown mesh for VMD')
 
@@ -109,11 +109,11 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         """Circular Loop, magnetic vector potential"""
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
-                MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
             else:
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
         elif mesh._meshType is 'TENSOR':
-            MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
         else:
             raise Exception('Unknown mesh for CircularLoop')
 

SimPEG/EM/Utils/__init__.py

@@ -1,5 +1,5 @@
 # import Sources
 # import Ana
 # import Solver
-import EMUtils
-import SrcUtils
+from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
+from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential

SimPEG/EM/Utils/testingUtils.py

@@ -0,0 +1,75 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+
+def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
+    cs = 5.
+    ncx, ncy, ncz = 6, 6, 6
+    npad = 3
+    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
+
+    mapping = Maps.ExpMap(mesh)
+
+    x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
+    XYZ = Utils.ndgrid(x,x,np.r_[0.])
+    Rx0 = EM.FDEM.Rx(XYZ, comp)
+
+    Src = []
+
+    for SrcType in SrcList:
+        if SrcType is 'MagDipole':
+            Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+        elif SrcType is 'MagDipole_Bfield':
+            Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+        elif SrcType is 'CircularLoop':
+            Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+        elif SrcType is 'RawVec':
+            if fdemType is 'e' or fdemType is 'b':
+                S_m = np.zeros(mesh.nF)
+                S_e = np.zeros(mesh.nE)
+                S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+                S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+
+            elif fdemType is 'h' or fdemType is 'j':
+                S_m = np.zeros(mesh.nE)
+                S_e = np.zeros(mesh.nF)
+                S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+                S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+
+    if verbose:
+        print '  Fetching %s problem' % (fdemType)
+
+    if fdemType == 'e':
+        survey = EM.FDEM.Survey(Src)
+        prb = EM.FDEM.Problem_e(mesh, mapping=mapping)
+
+    elif fdemType == 'b':
+        survey = EM.FDEM.Survey(Src)
+        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+
+    elif fdemType == 'j':
+        survey = EM.FDEM.Survey(Src)
+        prb = EM.FDEM.Problem_j(mesh, mapping=mapping)
+
+    elif fdemType == 'h':
+        survey = EM.FDEM.Survey(Src)
+        prb = EM.FDEM.Problem_h(mesh, mapping=mapping)
+
+    else:
+        raise NotImplementedError()
+    prb.pair(survey)
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        pass
+
+    return prb

SimPEG/Utils/__init__.py

@@ -7,4 +7,4 @@ from ipythonutils import easyAnimate as animate
 from CounterUtils import *
 import ModelBuilder
 import SolverUtils
-
+from coordutils import *

SimPEG/Utils/codeutils.py

@@ -3,10 +3,7 @@ import time
 import numpy as np
 from functools import wraps
 
-
-class SimPEGMetaClass(type):
-    def __new__(cls, name, bases, attrs):
-        return super(SimPEGMetaClass, cls).__new__(cls, name, bases, attrs)
+SimPEGMetaClass = type
 
 def memProfileWrapper(towrap, *funNames):
     """

SimPEG/Utils/coordutils.py

@@ -0,0 +1,62 @@
+import numpy as np
+from SimPEG.Utils import mkvc
+
+def rotationMatrixFromNormals(v0,v1,tol=1e-20):
+    """
+        Performs the minimum number of rotations to define a rotation from the direction indicated by the vector n0 to the direction indicated by n1.
+        The axis of rotation is n0 x n1
+        https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
+
+        :param numpy.array v0: vector of length 3
+        :param numpy.array v1: vector of length 3
+        :param tol = 1e-20: tolerance. If the norm of the cross product between the two vectors is below this, no rotation is performed
+        :rtype: numpy.array, 3x3
+        :return: rotation matrix which rotates the frame so that n0 is aligned with n1
+
+    """
+
+    # ensure both n0, n1 are vectors of length 1
+    assert len(v0) == 3, "Length of n0 should be 3"
+    assert len(v1) == 3, "Length of n1 should be 3"
+
+    # ensure both are true normals
+    n0 = v0*1./np.linalg.norm(v0)
+    n1 = v1*1./np.linalg.norm(v1)
+
+    n0dotn1 = n0.dot(n1) 
+
+    # define the rotation axis, which is the cross product of the two vectors
+    rotAx = np.cross(n0,n1)
+
+    if np.linalg.norm(rotAx) < tol:
+        return np.eye(3,dtype=float)
+
+    rotAx *= 1./np.linalg.norm(rotAx)
+
+    cosT = n0dotn1/(np.linalg.norm(n0)*np.linalg.norm(n1))
+    sinT = np.sqrt(1.-n0dotn1**2)
+
+    ux = np.array([[0., -rotAx[2], rotAx[1]], [rotAx[2], 0., -rotAx[0]], [-rotAx[1], rotAx[0], 0.]],dtype=float)
+
+    return np.eye(3,dtype=float) + sinT*ux + (1.-cosT)*(ux.dot(ux))
+
+
+def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]):
+    """
+        rotates a grid so that the vector n0 is aligned with the vector n1
+
+        :param numpy.array n0: vector of length 3, should have norm 1
+        :param numpy.array n1: vector of length 3, should have norm 1
+        :param numpy.array x0: vector of length 3, point about which we perform the rotation 
+        :rtype: numpy.array, 3x3
+        :return: rotation matrix which rotates the frame so that n0 is aligned with n1
+    """
+
+    R = rotationMatrixFromNormals(n0, n1)
+
+    assert XYZ.shape[1] == 3, "Grid XYZ should be 3 wide"
+    assert len(x0) == 3, "x0 should have length 3"
+
+    X0 = np.ones([XYZ.shape[0],1])*mkvc(x0)
+
+    return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0

tests/em/fdem/forward/__init__.py

@@ -0,0 +1,11 @@
+if __name__ == '__main__':
+    import os
+    import glob
+    import unittest
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/fdem/forward/test_FDEM_analytics.py

@@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase):
 
         x = np.linspace(-10,10,5)
         XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-        rxList = EM.FDEM.RxFDEM(XYZ, 'exi')
-        Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
+        rxList = EM.FDEM.Rx(XYZ, 'exi')
+        Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
 
-        survey = EM.FDEM.SurveyFDEM([Src0])
+        survey = EM.FDEM.Survey([Src0])
 
-        prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
         prb.pair(survey)
 
         try:
@@ -114,19 +114,19 @@ class FDEM_analyticTests(unittest.TestCase):
 
         de = np.zeros(mesh.nF,dtype=complex)
         de[s_ind] = 1./csz
-        de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)]
+        de_p = [EM.FDEM.Src.RawVec_e([],freq,de/mesh.area)]
 
-        dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)]
+        dm_p = [EM.FDEM.Src.MagDipole([],freq,src_loc)]
 
 
         # Pair the problem and survey
-        surveye = EM.FDEM.SurveyFDEM(de_p)
-        surveym = EM.FDEM.SurveyFDEM(dm_p)
+        surveye = EM.FDEM.Survey(de_p)
+        surveym = EM.FDEM.Survey(dm_p)
 
         mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
 
-        prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
-        prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
+        prbe = EM.FDEM.Problem_h(mesh, mapping=mapping)
+        prbm = EM.FDEM.Problem_e(mesh, mapping=mapping)
 
         prbe.pair(surveye) # pair problem and survey
         prbm.pair(surveym)

tests/em/fdem/forward/test_FDEM_forward.py

@@ -0,0 +1,127 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem
+
+testEB = True
+testHJ = True
+
+verbose = False
+
+TOL = 1e-5
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = 1e-1
+addrandoms = True
+
+SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
+
+
+def crossCheckTest(fdemType, comp):
+
+    l2norm = lambda r: np.sqrt(r.dot(r))
+
+    prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose)
+    mesh = prb1.mesh
+    print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
+    m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
+    mu = np.log(np.ones(mesh.nC)*MU)
+
+    if addrandoms is True:
+        m  = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
+        mu = mu + np.random.randn(mesh.nC)*MU*1e-1
+
+    # prb1.PropMap.PropModel.mu = mu
+    # prb1.PropMap.PropModel.mui = 1./mu
+    survey1 = prb1.survey
+    d1 = survey1.dpred(m)
+
+    if verbose:
+        print '  Problem 1 solved'
+
+    if fdemType == 'e':
+        prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose)
+    elif fdemType == 'b':
+        prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose)
+    elif fdemType == 'j':
+        prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose)
+    elif fdemType == 'h':
+        prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose)
+    else:
+        raise NotImplementedError()
+
+    # prb2.mu = mu
+    survey2 = prb2.survey
+    d2 = survey2.dpred(m)
+
+    if verbose:
+        print '  Problem 2 solved'
+
+    r = d2-d1
+    l2r = l2norm(r)
+
+    tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
+    print l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol
+    return l2r < tol
+
+
+class FDEM_CrossCheck(unittest.TestCase):
+    if testEB:
+        def test_EB_CrossCheck_exr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'exr'))
+        def test_EB_CrossCheck_eyr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'eyr'))
+        def test_EB_CrossCheck_ezr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'ezr'))
+        def test_EB_CrossCheck_exi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'exi'))
+        def test_EB_CrossCheck_eyi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'eyi'))
+        def test_EB_CrossCheck_ezi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'ezi'))
+
+        def test_EB_CrossCheck_bxr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'bxr'))
+        def test_EB_CrossCheck_byr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'byr'))
+        def test_EB_CrossCheck_bzr_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'bzr'))
+        def test_EB_CrossCheck_bxi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'bxi'))
+        def test_EB_CrossCheck_byi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'byi'))
+        def test_EB_CrossCheck_bzi_Eform(self):
+            self.assertTrue(crossCheckTest('e', 'bzi'))
+
+    if testHJ:
+        def test_HJ_CrossCheck_jxr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jxr'))
+        def test_HJ_CrossCheck_jyr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jyr'))
+        def test_HJ_CrossCheck_jzr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jzr'))
+        def test_HJ_CrossCheck_jxi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jxi'))
+        def test_HJ_CrossCheck_jyi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jyi'))
+        def test_HJ_CrossCheck_jzi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'jzi'))
+
+        def test_HJ_CrossCheck_hxr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hxr'))
+        def test_HJ_CrossCheck_hyr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hyr'))
+        def test_HJ_CrossCheck_hzr_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hzr'))
+        def test_HJ_CrossCheck_hxi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hxi'))
+        def test_HJ_CrossCheck_hyi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hyi'))
+        def test_HJ_CrossCheck_hzi_Jform(self):
+            self.assertTrue(crossCheckTest('j', 'hzi'))
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/fdem/inverse/adjoint/__init__.py

@@ -0,0 +1,11 @@
+if __name__ == '__main__':
+    import os
+    import glob
+    import unittest
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py

@@ -0,0 +1,156 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem
+
+testEB = True
+testHJ = True
+
+verbose = False
+
+TOL = 1e-5
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = 1e-1
+addrandoms = True
+
+SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+
+def adjointTest(fdemType, comp):
+    prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
+    print 'Adjoint %s formulation - %s' % (fdemType, comp)
+
+    m  = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
+    mu = np.ones(prb.mesh.nC)*MU
+
+    if addrandoms is True:
+        m  = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
+        mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
+
+    survey = prb.survey
+    # prb.PropMap.PropModel.mu = mu
+    # prb.PropMap.PropModel.mui = 1./mu
+    u = prb.fields(m)
+
+    v = np.random.rand(survey.nD)
+    w = np.random.rand(prb.mesh.nC)
+
+    vJw = v.dot(prb.Jvec(m, w, u))
+    wJtv = w.dot(prb.Jtvec(m, v, u))
+    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
+    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
+    return np.abs(vJw - wJtv) < tol
+
+class FDEM_AdjointTests(unittest.TestCase):
+    if testEB:
+        def test_Jtvec_adjointTest_exr_Eform(self):
+            self.assertTrue(adjointTest('e', 'exr'))
+        def test_Jtvec_adjointTest_eyr_Eform(self):
+            self.assertTrue(adjointTest('e', 'eyr'))
+        def test_Jtvec_adjointTest_ezr_Eform(self):
+            self.assertTrue(adjointTest('e', 'ezr'))
+        def test_Jtvec_adjointTest_exi_Eform(self):
+            self.assertTrue(adjointTest('e', 'exi'))
+        def test_Jtvec_adjointTest_eyi_Eform(self):
+            self.assertTrue(adjointTest('e', 'eyi'))
+        def test_Jtvec_adjointTest_ezi_Eform(self):
+            self.assertTrue(adjointTest('e', 'ezi'))
+
+        def test_Jtvec_adjointTest_bxr_Eform(self):
+            self.assertTrue(adjointTest('e', 'bxr'))
+        def test_Jtvec_adjointTest_byr_Eform(self):
+            self.assertTrue(adjointTest('e', 'byr'))
+        def test_Jtvec_adjointTest_bzr_Eform(self):
+            self.assertTrue(adjointTest('e', 'bzr'))
+        def test_Jtvec_adjointTest_bxi_Eform(self):
+            self.assertTrue(adjointTest('e', 'bxi'))
+        def test_Jtvec_adjointTest_byi_Eform(self):
+            self.assertTrue(adjointTest('e', 'byi'))
+        def test_Jtvec_adjointTest_bzi_Eform(self):
+            self.assertTrue(adjointTest('e', 'bzi'))
+
+        def test_Jtvec_adjointTest_exr_Bform(self):
+            self.assertTrue(adjointTest('b', 'exr'))
+        def test_Jtvec_adjointTest_eyr_Bform(self):
+            self.assertTrue(adjointTest('b', 'eyr'))
+        def test_Jtvec_adjointTest_ezr_Bform(self):
+            self.assertTrue(adjointTest('b', 'ezr'))
+        def test_Jtvec_adjointTest_exi_Bform(self):
+            self.assertTrue(adjointTest('b', 'exi'))
+        def test_Jtvec_adjointTest_eyi_Bform(self):
+            self.assertTrue(adjointTest('b', 'eyi'))
+        def test_Jtvec_adjointTest_ezi_Bform(self):
+            self.assertTrue(adjointTest('b', 'ezi'))
+        def test_Jtvec_adjointTest_bxr_Bform(self):
+            self.assertTrue(adjointTest('b', 'bxr'))
+        def test_Jtvec_adjointTest_byr_Bform(self):
+            self.assertTrue(adjointTest('b', 'byr'))
+        def test_Jtvec_adjointTest_bzr_Bform(self):
+            self.assertTrue(adjointTest('b', 'bzr'))
+        def test_Jtvec_adjointTest_bxi_Bform(self):
+            self.assertTrue(adjointTest('b', 'bxi'))
+        def test_Jtvec_adjointTest_byi_Bform(self):
+            self.assertTrue(adjointTest('b', 'byi'))
+        def test_Jtvec_adjointTest_bzi_Bform(self):
+            self.assertTrue(adjointTest('b', 'bzi'))
+
+
+    if testHJ:
+        def test_Jtvec_adjointTest_jxr_Jform(self):
+            self.assertTrue(adjointTest('j', 'jxr'))
+        def test_Jtvec_adjointTest_jyr_Jform(self):
+            self.assertTrue(adjointTest('j', 'jyr'))
+        def test_Jtvec_adjointTest_jzr_Jform(self):
+            self.assertTrue(adjointTest('j', 'jzr'))
+        def test_Jtvec_adjointTest_jxi_Jform(self):
+            self.assertTrue(adjointTest('j', 'jxi'))
+        def test_Jtvec_adjointTest_jyi_Jform(self):
+            self.assertTrue(adjointTest('j', 'jyi'))
+        def test_Jtvec_adjointTest_jzi_Jform(self):
+            self.assertTrue(adjointTest('j', 'jzi'))
+
+        def test_Jtvec_adjointTest_hxr_Jform(self):
+            self.assertTrue(adjointTest('j', 'hxr'))
+        def test_Jtvec_adjointTest_hyr_Jform(self):
+            self.assertTrue(adjointTest('j', 'hyr'))
+        def test_Jtvec_adjointTest_hzr_Jform(self):
+            self.assertTrue(adjointTest('j', 'hzr'))
+        def test_Jtvec_adjointTest_hxi_Jform(self):
+            self.assertTrue(adjointTest('j', 'hxi'))
+        def test_Jtvec_adjointTest_hyi_Jform(self):
+            self.assertTrue(adjointTest('j', 'hyi'))
+        def test_Jtvec_adjointTest_hzi_Jform(self):
+             self.assertTrue(adjointTest('j', 'hzi'))
+
+        def test_Jtvec_adjointTest_hxr_Hform(self):
+            self.assertTrue(adjointTest('h', 'hxr'))
+        def test_Jtvec_adjointTest_hyr_Hform(self):
+            self.assertTrue(adjointTest('h', 'hyr'))
+        def test_Jtvec_adjointTest_hzr_Hform(self):
+            self.assertTrue(adjointTest('h', 'hzr'))
+        def test_Jtvec_adjointTest_hxi_Hform(self):
+            self.assertTrue(adjointTest('h', 'hxi'))
+        def test_Jtvec_adjointTest_hyi_Hform(self):
+            self.assertTrue(adjointTest('h', 'hyi'))
+        def test_Jtvec_adjointTest_hzi_Hform(self):
+            self.assertTrue(adjointTest('h', 'hzi'))
+
+        def test_Jtvec_adjointTest_hxr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jxr'))
+        def test_Jtvec_adjointTest_hyr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jyr'))
+        def test_Jtvec_adjointTest_hzr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jzr'))
+        def test_Jtvec_adjointTest_hxi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jxi'))
+        def test_Jtvec_adjointTest_hyi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jyi'))
+        def test_Jtvec_adjointTest_hzi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jzi'))
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/fdem/inverse/derivs/__init__.py

@@ -0,0 +1,11 @@
+if __name__ == '__main__':
+    import os
+    import glob
+    import unittest
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/fdem/inverse/derivs/test_FDEM_derivs.py

@@ -0,0 +1,154 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem
+
+testDerivs = True
+testEB = True
+testHJ = True
+
+verbose = False
+
+TOL = 1e-5
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = 1e-1
+addrandoms = True
+
+SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+
+
+def derivTest(fdemType, comp):
+
+    prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
+    print '%s formulation - %s' % (fdemType, comp)
+    x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
+    mu = np.log(np.ones(prb.mesh.nC)*MU)
+
+    if addrandoms is True:
+        x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
+        mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
+
+    # prb.PropMap.PropModel.mu = mu
+    # prb.PropMap.PropModel.mui = 1./mu
+
+    survey = prb.survey
+    def fun(x):
+        return survey.dpred(x), lambda x: prb.Jvec(x0, x)
+    return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
+
+
+class FDEM_DerivTests(unittest.TestCase):
+
+    if testEB:
+        def test_Jvec_exr_Eform(self):
+            self.assertTrue(derivTest('e', 'exr'))
+        def test_Jvec_eyr_Eform(self):
+            self.assertTrue(derivTest('e', 'eyr'))
+        def test_Jvec_ezr_Eform(self):
+            self.assertTrue(derivTest('e', 'ezr'))
+        def test_Jvec_exi_Eform(self):
+            self.assertTrue(derivTest('e', 'exi'))
+        def test_Jvec_eyi_Eform(self):
+            self.assertTrue(derivTest('e', 'eyi'))
+        def test_Jvec_ezi_Eform(self):
+            self.assertTrue(derivTest('e', 'ezi'))
+
+        def test_Jvec_bxr_Eform(self):
+            self.assertTrue(derivTest('e', 'bxr'))
+        def test_Jvec_byr_Eform(self):
+            self.assertTrue(derivTest('e', 'byr'))
+        def test_Jvec_bzr_Eform(self):
+            self.assertTrue(derivTest('e', 'bzr'))
+        def test_Jvec_bxi_Eform(self):
+            self.assertTrue(derivTest('e', 'bxi'))
+        def test_Jvec_byi_Eform(self):
+            self.assertTrue(derivTest('e', 'byi'))
+        def test_Jvec_bzi_Eform(self):
+            self.assertTrue(derivTest('e', 'bzi'))
+
+        def test_Jvec_exr_Bform(self):
+            self.assertTrue(derivTest('b', 'exr'))
+        def test_Jvec_eyr_Bform(self):
+            self.assertTrue(derivTest('b', 'eyr'))
+        def test_Jvec_ezr_Bform(self):
+            self.assertTrue(derivTest('b', 'ezr'))
+        def test_Jvec_exi_Bform(self):
+            self.assertTrue(derivTest('b', 'exi'))
+        def test_Jvec_eyi_Bform(self):
+            self.assertTrue(derivTest('b', 'eyi'))
+        def test_Jvec_ezi_Bform(self):
+            self.assertTrue(derivTest('b', 'ezi'))
+
+        def test_Jvec_bxr_Bform(self):
+            self.assertTrue(derivTest('b', 'bxr'))
+        def test_Jvec_byr_Bform(self):
+            self.assertTrue(derivTest('b', 'byr'))
+        def test_Jvec_bzr_Bform(self):
+            self.assertTrue(derivTest('b', 'bzr'))
+        def test_Jvec_bxi_Bform(self):
+            self.assertTrue(derivTest('b', 'bxi'))
+        def test_Jvec_byi_Bform(self):
+            self.assertTrue(derivTest('b', 'byi'))
+        def test_Jvec_bzi_Bform(self):
+            self.assertTrue(derivTest('b', 'bzi'))
+
+    if testHJ:
+        def test_Jvec_jxr_Jform(self):
+            self.assertTrue(derivTest('j', 'jxr'))
+        def test_Jvec_jyr_Jform(self):
+            self.assertTrue(derivTest('j', 'jyr'))
+        def test_Jvec_jzr_Jform(self):
+            self.assertTrue(derivTest('j', 'jzr'))
+        def test_Jvec_jxi_Jform(self):
+            self.assertTrue(derivTest('j', 'jxi'))
+        def test_Jvec_jyi_Jform(self):
+            self.assertTrue(derivTest('j', 'jyi'))
+        def test_Jvec_jzi_Jform(self):
+            self.assertTrue(derivTest('j', 'jzi'))
+
+        def test_Jvec_hxr_Jform(self):
+            self.assertTrue(derivTest('j', 'hxr'))
+        def test_Jvec_hyr_Jform(self):
+            self.assertTrue(derivTest('j', 'hyr'))
+        def test_Jvec_hzr_Jform(self):
+            self.assertTrue(derivTest('j', 'hzr'))
+        def test_Jvec_hxi_Jform(self):
+            self.assertTrue(derivTest('j', 'hxi'))
+        def test_Jvec_hyi_Jform(self):
+            self.assertTrue(derivTest('j', 'hyi'))
+        def test_Jvec_hzi_Jform(self):
+            self.assertTrue(derivTest('j', 'hzi'))
+
+        def test_Jvec_hxr_Hform(self):
+            self.assertTrue(derivTest('h', 'hxr'))
+        def test_Jvec_hyr_Hform(self):
+            self.assertTrue(derivTest('h', 'hyr'))
+        def test_Jvec_hzr_Hform(self):
+            self.assertTrue(derivTest('h', 'hzr'))
+        def test_Jvec_hxi_Hform(self):
+            self.assertTrue(derivTest('h', 'hxi'))
+        def test_Jvec_hyi_Hform(self):
+            self.assertTrue(derivTest('h', 'hyi'))
+        def test_Jvec_hzi_Hform(self):
+            self.assertTrue(derivTest('h', 'hzi'))
+
+        def test_Jvec_hxr_Hform(self):
+            self.assertTrue(derivTest('h', 'jxr'))
+        def test_Jvec_hyr_Hform(self):
+            self.assertTrue(derivTest('h', 'jyr'))
+        def test_Jvec_hzr_Hform(self):
+            self.assertTrue(derivTest('h', 'jzr'))
+        def test_Jvec_hxi_Hform(self):
+            self.assertTrue(derivTest('h', 'jxi'))
+        def test_Jvec_hyi_Hform(self):
+            self.assertTrue(derivTest('h', 'jyi'))
+        def test_Jvec_hzi_Hform(self):
+            self.assertTrue(derivTest('h', 'jzi'))
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/tdem/__init__.py

@@ -0,0 +1,11 @@
+if __name__ == '__main__':
+    import os
+    import glob
+    import unittest
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/test_FDEM.py

@@ -1,478 +0,0 @@
-import unittest
-from SimPEG import *
-from SimPEG import EM
-import sys
-from scipy.constants import mu_0
-
-testDerivs = True
-testCrossCheck = True
-testAdjoint = True
-testEB = True
-testHJ = True
-
-verbose = False
-
-TOL = 1e-5
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = 1e-1
-addrandoms = True
-
-SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
-
-
-def getProblem(fdemType, comp):
-    cs = 5.
-    ncx, ncy, ncz = 6, 6, 6
-    npad = 3
-    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
-    mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
-
-    mapping = Maps.ExpMap(mesh)
-
-    x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
-    XYZ = Utils.ndgrid(x,x,np.r_[0.])
-    Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
-
-    if SrcType is 'MagDipole':
-        Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
-    elif SrcType is 'MagDipole_Bfield':
-        Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
-    elif SrcType is 'CircularLoop':
-        Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
-
-    if verbose:
-        print '  Fetching %s problem' % (fdemType)
-
-    if fdemType == 'e':
-        if SrcType is 'RawVec':
-            S_m = np.zeros(mesh.nF)
-            S_e = np.zeros(mesh.nE)
-            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
-            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
-            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
-
-        survey = EM.FDEM.SurveyFDEM([Src])
-        prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
-
-    elif fdemType == 'b':
-        if SrcType is 'RawVec':
-            S_m = np.zeros(mesh.nF)
-            S_e = np.zeros(mesh.nE)
-            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
-            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
-            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
-
-        survey = EM.FDEM.SurveyFDEM([Src])
-        prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
-
-    elif fdemType == 'j':
-        if SrcType is 'RawVec':
-            S_m = np.zeros(mesh.nE)
-            S_e = np.zeros(mesh.nF)
-            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
-            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
-            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
-
-        survey = EM.FDEM.SurveyFDEM([Src])
-        prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
-
-    elif fdemType == 'h':
-        if SrcType is 'RawVec':
-            S_m = np.zeros(mesh.nE)
-            S_e = np.zeros(mesh.nF)
-            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
-            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
-            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
-
-        survey = EM.FDEM.SurveyFDEM([Src])
-        prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
-
-    else:
-        raise NotImplementedError()
-    prb.pair(survey)
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        pass
-
-    return prb
-
-def adjointTest(fdemType, comp):
-    prb = getProblem(fdemType, comp)
-    print 'Adjoint %s formulation - %s' % (fdemType, comp)
-
-    m  = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
-    mu = np.ones(prb.mesh.nC)*MU
-
-    if addrandoms is True:
-        m  = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
-        mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
-
-    survey = prb.survey
-    # prb.PropMap.PropModel.mu = mu
-    # prb.PropMap.PropModel.mui = 1./mu
-    u = prb.fields(m)
-
-    v = np.random.rand(survey.nD)
-    w = np.random.rand(prb.mesh.nC)
-
-    vJw = v.dot(prb.Jvec(m, w, u))
-    wJtv = w.dot(prb.Jtvec(m, v, u))
-    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
-    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
-    return np.abs(vJw - wJtv) < tol
-
-
-def derivTest(fdemType, comp):
-
-    prb = getProblem(fdemType, comp)
-    print '%s formulation - %s' % (fdemType, comp)
-    x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
-    mu = np.log(np.ones(prb.mesh.nC)*MU)
-
-    if addrandoms is True:
-        x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
-        mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
-
-    # prb.PropMap.PropModel.mu = mu
-    # prb.PropMap.PropModel.mui = 1./mu
-
-    survey = prb.survey
-    def fun(x):
-        return survey.dpred(x), lambda x: prb.Jvec(x0, x)
-    return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
-
-
-def crossCheckTest(fdemType, comp):
-
-    l2norm = lambda r: np.sqrt(r.dot(r))
-
-    prb1 = getProblem(fdemType, comp)
-    mesh = prb1.mesh
-    print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
-    m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
-    mu = np.log(np.ones(mesh.nC)*MU)
-
-    if addrandoms is True:
-        m  = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
-        mu = mu + np.random.randn(mesh.nC)*MU*1e-1
-
-    # prb1.PropMap.PropModel.mu = mu
-    # prb1.PropMap.PropModel.mui = 1./mu
-    survey1 = prb1.survey
-    d1 = survey1.dpred(m)
-
-    if verbose:
-        print '  Problem 1 solved'
-
-    if fdemType == 'e':
-        prb2 = getProblem('b', comp)
-    elif fdemType == 'b':
-        prb2 = getProblem('e', comp)
-    elif fdemType == 'j':
-        prb2 = getProblem('h', comp)
-    elif fdemType == 'h':
-        prb2 = getProblem('j', comp)
-    else:
-        raise NotImplementedError()
-
-    # prb2.mu = mu
-    survey2 = prb2.survey
-    d2 = survey2.dpred(m)
-
-    if verbose:
-        print '  Problem 2 solved'
-
-    r = d2-d1
-    l2r = l2norm(r)
-
-    tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
-    print l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol
-    return l2r < tol
-
-
-
-
-
-class FDEM_DerivTests(unittest.TestCase):
-
-    if testDerivs:
-        if testEB:
-            def test_Jvec_exr_Eform(self):
-                self.assertTrue(derivTest('e', 'exr'))
-            def test_Jvec_eyr_Eform(self):
-                self.assertTrue(derivTest('e', 'eyr'))
-            def test_Jvec_ezr_Eform(self):
-                self.assertTrue(derivTest('e', 'ezr'))
-            def test_Jvec_exi_Eform(self):
-                self.assertTrue(derivTest('e', 'exi'))
-            def test_Jvec_eyi_Eform(self):
-                self.assertTrue(derivTest('e', 'eyi'))
-            def test_Jvec_ezi_Eform(self):
-                self.assertTrue(derivTest('e', 'ezi'))
-
-            def test_Jvec_bxr_Eform(self):
-                self.assertTrue(derivTest('e', 'bxr'))
-            def test_Jvec_byr_Eform(self):
-                self.assertTrue(derivTest('e', 'byr'))
-            def test_Jvec_bzr_Eform(self):
-                self.assertTrue(derivTest('e', 'bzr'))
-            def test_Jvec_bxi_Eform(self):
-                self.assertTrue(derivTest('e', 'bxi'))
-            def test_Jvec_byi_Eform(self):
-                self.assertTrue(derivTest('e', 'byi'))
-            def test_Jvec_bzi_Eform(self):
-                self.assertTrue(derivTest('e', 'bzi'))
-
-            def test_Jvec_exr_Bform(self):
-                self.assertTrue(derivTest('b', 'exr'))
-            def test_Jvec_eyr_Bform(self):
-                self.assertTrue(derivTest('b', 'eyr'))
-            def test_Jvec_ezr_Bform(self):
-                self.assertTrue(derivTest('b', 'ezr'))
-            def test_Jvec_exi_Bform(self):
-                self.assertTrue(derivTest('b', 'exi'))
-            def test_Jvec_eyi_Bform(self):
-                self.assertTrue(derivTest('b', 'eyi'))
-            def test_Jvec_ezi_Bform(self):
-                self.assertTrue(derivTest('b', 'ezi'))
-
-            def test_Jvec_bxr_Bform(self):
-                self.assertTrue(derivTest('b', 'bxr'))
-            def test_Jvec_byr_Bform(self):
-                self.assertTrue(derivTest('b', 'byr'))
-            def test_Jvec_bzr_Bform(self):
-                self.assertTrue(derivTest('b', 'bzr'))
-            def test_Jvec_bxi_Bform(self):
-                self.assertTrue(derivTest('b', 'bxi'))
-            def test_Jvec_byi_Bform(self):
-                self.assertTrue(derivTest('b', 'byi'))
-            def test_Jvec_bzi_Bform(self):
-                self.assertTrue(derivTest('b', 'bzi'))
-
-        if testHJ:
-            def test_Jvec_jxr_Jform(self):
-                self.assertTrue(derivTest('j', 'jxr'))
-            def test_Jvec_jyr_Jform(self):
-                self.assertTrue(derivTest('j', 'jyr'))
-            def test_Jvec_jzr_Jform(self):
-                self.assertTrue(derivTest('j', 'jzr'))
-            def test_Jvec_jxi_Jform(self):
-                self.assertTrue(derivTest('j', 'jxi'))
-            def test_Jvec_jyi_Jform(self):
-                self.assertTrue(derivTest('j', 'jyi'))
-            def test_Jvec_jzi_Jform(self):
-                self.assertTrue(derivTest('j', 'jzi'))
-
-            def test_Jvec_hxr_Jform(self):
-                self.assertTrue(derivTest('j', 'hxr'))
-            def test_Jvec_hyr_Jform(self):
-                self.assertTrue(derivTest('j', 'hyr'))
-            def test_Jvec_hzr_Jform(self):
-                self.assertTrue(derivTest('j', 'hzr'))
-            def test_Jvec_hxi_Jform(self):
-                self.assertTrue(derivTest('j', 'hxi'))
-            def test_Jvec_hyi_Jform(self):
-                self.assertTrue(derivTest('j', 'hyi'))
-            def test_Jvec_hzi_Jform(self):
-                self.assertTrue(derivTest('j', 'hzi'))
-
-            def test_Jvec_hxr_Hform(self):
-                self.assertTrue(derivTest('h', 'hxr'))
-            def test_Jvec_hyr_Hform(self):
-                self.assertTrue(derivTest('h', 'hyr'))
-            def test_Jvec_hzr_Hform(self):
-                self.assertTrue(derivTest('h', 'hzr'))
-            def test_Jvec_hxi_Hform(self):
-                self.assertTrue(derivTest('h', 'hxi'))
-            def test_Jvec_hyi_Hform(self):
-                self.assertTrue(derivTest('h', 'hyi'))
-            def test_Jvec_hzi_Hform(self):
-                self.assertTrue(derivTest('h', 'hzi'))
-
-            def test_Jvec_hxr_Hform(self):
-                self.assertTrue(derivTest('h', 'jxr'))
-            def test_Jvec_hyr_Hform(self):
-                self.assertTrue(derivTest('h', 'jyr'))
-            def test_Jvec_hzr_Hform(self):
-                self.assertTrue(derivTest('h', 'jzr'))
-            def test_Jvec_hxi_Hform(self):
-                self.assertTrue(derivTest('h', 'jxi'))
-            def test_Jvec_hyi_Hform(self):
-                self.assertTrue(derivTest('h', 'jyi'))
-            def test_Jvec_hzi_Hform(self):
-                self.assertTrue(derivTest('h', 'jzi'))
-
-
-    if testAdjoint:
-        if testEB:
-            def test_Jtvec_adjointTest_exr_Eform(self):
-                self.assertTrue(adjointTest('e', 'exr'))
-            def test_Jtvec_adjointTest_eyr_Eform(self):
-                self.assertTrue(adjointTest('e', 'eyr'))
-            def test_Jtvec_adjointTest_ezr_Eform(self):
-                self.assertTrue(adjointTest('e', 'ezr'))
-            def test_Jtvec_adjointTest_exi_Eform(self):
-                self.assertTrue(adjointTest('e', 'exi'))
-            def test_Jtvec_adjointTest_eyi_Eform(self):
-                self.assertTrue(adjointTest('e', 'eyi'))
-            def test_Jtvec_adjointTest_ezi_Eform(self):
-                self.assertTrue(adjointTest('e', 'ezi'))
-
-            def test_Jtvec_adjointTest_bxr_Eform(self):
-                self.assertTrue(adjointTest('e', 'bxr'))
-            def test_Jtvec_adjointTest_byr_Eform(self):
-                self.assertTrue(adjointTest('e', 'byr'))
-            def test_Jtvec_adjointTest_bzr_Eform(self):
-                self.assertTrue(adjointTest('e', 'bzr'))
-            def test_Jtvec_adjointTest_bxi_Eform(self):
-                self.assertTrue(adjointTest('e', 'bxi'))
-            def test_Jtvec_adjointTest_byi_Eform(self):
-                self.assertTrue(adjointTest('e', 'byi'))
-            def test_Jtvec_adjointTest_bzi_Eform(self):
-                self.assertTrue(adjointTest('e', 'bzi'))
-
-            def test_Jtvec_adjointTest_exr_Bform(self):
-                self.assertTrue(adjointTest('b', 'exr'))
-            def test_Jtvec_adjointTest_eyr_Bform(self):
-                self.assertTrue(adjointTest('b', 'eyr'))
-            def test_Jtvec_adjointTest_ezr_Bform(self):
-                self.assertTrue(adjointTest('b', 'ezr'))
-            def test_Jtvec_adjointTest_exi_Bform(self):
-                self.assertTrue(adjointTest('b', 'exi'))
-            def test_Jtvec_adjointTest_eyi_Bform(self):
-                self.assertTrue(adjointTest('b', 'eyi'))
-            def test_Jtvec_adjointTest_ezi_Bform(self):
-                self.assertTrue(adjointTest('b', 'ezi'))
-            def test_Jtvec_adjointTest_bxr_Bform(self):
-                self.assertTrue(adjointTest('b', 'bxr'))
-            def test_Jtvec_adjointTest_byr_Bform(self):
-                self.assertTrue(adjointTest('b', 'byr'))
-            def test_Jtvec_adjointTest_bzr_Bform(self):
-                self.assertTrue(adjointTest('b', 'bzr'))
-            def test_Jtvec_adjointTest_bxi_Bform(self):
-                self.assertTrue(adjointTest('b', 'bxi'))
-            def test_Jtvec_adjointTest_byi_Bform(self):
-                self.assertTrue(adjointTest('b', 'byi'))
-            def test_Jtvec_adjointTest_bzi_Bform(self):
-                self.assertTrue(adjointTest('b', 'bzi'))
-
-
-        if testHJ:
-            def test_Jtvec_adjointTest_jxr_Jform(self):
-                self.assertTrue(adjointTest('j', 'jxr'))
-            def test_Jtvec_adjointTest_jyr_Jform(self):
-                self.assertTrue(adjointTest('j', 'jyr'))
-            def test_Jtvec_adjointTest_jzr_Jform(self):
-                self.assertTrue(adjointTest('j', 'jzr'))
-            def test_Jtvec_adjointTest_jxi_Jform(self):
-                self.assertTrue(adjointTest('j', 'jxi'))
-            def test_Jtvec_adjointTest_jyi_Jform(self):
-                self.assertTrue(adjointTest('j', 'jyi'))
-            def test_Jtvec_adjointTest_jzi_Jform(self):
-                self.assertTrue(adjointTest('j', 'jzi'))
-
-            def test_Jtvec_adjointTest_hxr_Jform(self):
-                self.assertTrue(adjointTest('j', 'hxr'))
-            def test_Jtvec_adjointTest_hyr_Jform(self):
-                self.assertTrue(adjointTest('j', 'hyr'))
-            def test_Jtvec_adjointTest_hzr_Jform(self):
-                self.assertTrue(adjointTest('j', 'hzr'))
-            def test_Jtvec_adjointTest_hxi_Jform(self):
-                self.assertTrue(adjointTest('j', 'hxi'))
-            def test_Jtvec_adjointTest_hyi_Jform(self):
-                self.assertTrue(adjointTest('j', 'hyi'))
-            def test_Jtvec_adjointTest_hzi_Jform(self):
-                 self.assertTrue(adjointTest('j', 'hzi'))
-
-            def test_Jtvec_adjointTest_hxr_Hform(self):
-                self.assertTrue(adjointTest('h', 'hxr'))
-            def test_Jtvec_adjointTest_hyr_Hform(self):
-                self.assertTrue(adjointTest('h', 'hyr'))
-            def test_Jtvec_adjointTest_hzr_Hform(self):
-                self.assertTrue(adjointTest('h', 'hzr'))
-            def test_Jtvec_adjointTest_hxi_Hform(self):
-                self.assertTrue(adjointTest('h', 'hxi'))
-            def test_Jtvec_adjointTest_hyi_Hform(self):
-                self.assertTrue(adjointTest('h', 'hyi'))
-            def test_Jtvec_adjointTest_hzi_Hform(self):
-                self.assertTrue(adjointTest('h', 'hzi'))
-
-            def test_Jtvec_adjointTest_hxr_Hform(self):
-                self.assertTrue(adjointTest('h', 'jxr'))
-            def test_Jtvec_adjointTest_hyr_Hform(self):
-                self.assertTrue(adjointTest('h', 'jyr'))
-            def test_Jtvec_adjointTest_hzr_Hform(self):
-                self.assertTrue(adjointTest('h', 'jzr'))
-            def test_Jtvec_adjointTest_hxi_Hform(self):
-                self.assertTrue(adjointTest('h', 'jxi'))
-            def test_Jtvec_adjointTest_hyi_Hform(self):
-                self.assertTrue(adjointTest('h', 'jyi'))
-            def test_Jtvec_adjointTest_hzi_Hform(self):
-                self.assertTrue(adjointTest('h', 'jzi'))
-
-
-    if testCrossCheck:
-        if testEB:
-            def test_EB_CrossCheck_exr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'exr'))
-            def test_EB_CrossCheck_eyr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'eyr'))
-            def test_EB_CrossCheck_ezr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'ezr'))
-            def test_EB_CrossCheck_exi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'exi'))
-            def test_EB_CrossCheck_eyi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'eyi'))
-            def test_EB_CrossCheck_ezi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'ezi'))
-
-            def test_EB_CrossCheck_bxr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'bxr'))
-            def test_EB_CrossCheck_byr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'byr'))
-            def test_EB_CrossCheck_bzr_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'bzr'))
-            def test_EB_CrossCheck_bxi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'bxi'))
-            def test_EB_CrossCheck_byi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'byi'))
-            def test_EB_CrossCheck_bzi_Eform(self):
-                self.assertTrue(crossCheckTest('e', 'bzi'))
-
-        if testHJ:
-            def test_HJ_CrossCheck_jxr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jxr'))
-            def test_HJ_CrossCheck_jyr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jyr'))
-            def test_HJ_CrossCheck_jzr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jzr'))
-            def test_HJ_CrossCheck_jxi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jxi'))
-            def test_HJ_CrossCheck_jyi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jyi'))
-            def test_HJ_CrossCheck_jzi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'jzi'))
-
-            def test_HJ_CrossCheck_hxr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hxr'))
-            def test_HJ_CrossCheck_hyr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hyr'))
-            def test_HJ_CrossCheck_hzr_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hzr'))
-            def test_HJ_CrossCheck_hxi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hxi'))
-            def test_HJ_CrossCheck_hyi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hyi'))
-            def test_HJ_CrossCheck_hzi_Jform(self):
-                self.assertTrue(crossCheckTest('j', 'hzi'))
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/utils/test_coordutils.py

@@ -0,0 +1,56 @@
+import unittest, os
+import numpy as np
+from SimPEG import Utils
+
+tol = 1e-15
+
+class coorUtilsTest(unittest.TestCase):
+
+    def test_rotationMatrixFromNormals(self):
+        np.random.seed(0)
+        v0 = np.random.rand(3)
+        v0 *= 1./np.linalg.norm(v0)
+
+        np.random.seed(5)
+        v1 = np.random.rand(3)
+        v1 *= 1./np.linalg.norm(v1)
+
+        Rf = Utils.coordutils.rotationMatrixFromNormals(v0,v1)
+        Ri = Utils.coordutils.rotationMatrixFromNormals(v1,v0)
+
+        self.assertTrue(np.linalg.norm(Utils.mkvc(Rf.dot(v0) - v1)) < tol)
+        self.assertTrue(np.linalg.norm(Utils.mkvc(Ri.dot(v1) - v0)) < tol)
+
+    def test_rotatePointsFromNormals(self):
+        np.random.seed(10)
+        v0 = np.random.rand(3)
+        v0*= 1./np.linalg.norm(v0)
+
+        np.random.seed(15)
+        v1 = np.random.rand(3) 
+        v1*= 1./np.linalg.norm(v1)   
+
+        v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1)) 
+        
+        self.assertTrue(np.linalg.norm(v2-v1) < tol)
+
+    def test_rotateMatrixFromNormals(self):
+        np.random.seed(20)
+        n0 = np.random.rand(3)
+        n0 *= 1./np.linalg.norm(n0)
+
+        np.random.seed(25)
+        n1 = np.random.rand(3) 
+        n1 *= 1./np.linalg.norm(n1) 
+
+        np.random.seed(30)
+        scale =  np.random.rand(100,1)
+        XYZ0 = scale * n0
+        XYZ1 = scale * n1
+
+        XYZ2 = Utils.coordutils.rotatePointsFromNormals(XYZ0,n0,n1)
+        self.assertTrue(np.linalg.norm(Utils.mkvc(XYZ1) - Utils.mkvc(XYZ2))/np.linalg.norm(Utils.mkvc(XYZ1)) < tol)
+
+if __name__ == '__main__':
+    unittest.main()
+