diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py
index 603d322c..e4e30dcf 100644
--- a/SimPEG/EM/FDEM/FDEM.py
+++ b/SimPEG/EM/FDEM/FDEM.py
@@ -138,14 +138,14 @@ class BaseFDEMProblem(BaseEMProblem):
                     dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True)
                     du_dmT = -dA_dmT + dRHS_dmT
 
-                    df_dmT += du_dmT
+                    Df_DmT = df_dmT + du_dmT
 
                     # TODO: this should be taken care of by the reciever?
                     real_or_imag = rx.projComp
                     if real_or_imag is 'real':
-                        Jtv +=   np.array(df_dmT,dtype=complex).real
+                        Jtv +=   np.array(Df_DmT,dtype=complex).real
                     elif real_or_imag is 'imag':
-                        Jtv += - np.array(df_dmT,dtype=complex).real
+                        Jtv += - np.array(Df_DmT,dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
             
diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index 8d8c16f4..37a13f24 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -3,7 +3,7 @@ import scipy.sparse as sp
 import SimPEG
 from SimPEG import Utils
 from SimPEG.EM.Utils import omega
-from SimPEG.Utils import Zero, Identity
+from SimPEG.Utils import Zero, Identity, sdiag
 
 
 class Fields(SimPEG.Problem.Fields):
@@ -176,15 +176,34 @@ class Fields_e(Fields):
                     'eSecondary' : ['eSolution','E','_eSecondary'],
                     'b' : ['eSolution','F','_b'],
                     'bPrimary' : ['eSolution','F','_bPrimary'],
-                    'bSecondary' : ['eSolution','F','_bSecondary']
+                    'bSecondary' : ['eSolution','F','_bSecondary'],
+                    'j' : ['eSolution','CCV','_j'],
+                    'h' : ['eSolution','CCV','_h'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
+    def __init__(self, mesh, survey, **kwargs):
         Fields.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
+        self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
+        self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
+        self._nC = self.survey.prob.mesh.nC
+        self._MeSigma = self.survey.prob.MeSigma
+        self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
+        self._MfMui = self.survey.prob.MfMui
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
 
     def _ePrimary(self, eSolution, srcList):
         """
@@ -196,7 +215,7 @@ class Fields_e(Fields):
         :return: primary electric field as defined by the sources
         """
 
-        ePrimary = np.zeros_like(eSolution)
+        ePrimary = np.zeros([self.prob.mesh.nE,len(srcList)], dtype = complex)
         for i, src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
             ePrimary[:,i] = ePrimary[:,i] + ep
@@ -211,10 +230,8 @@ class Fields_e(Fields):
         :rtype: numpy.ndarray
         :return: secondary electric field
         """
-
         return eSolution
 
-
     def _eDeriv_u(self, src, v, adjoint = False):
         """
         Partial derivative of the total electric field with respect to the thing we 
@@ -253,7 +270,7 @@ class Fields_e(Fields):
         :return: primary magnetic flux density as defined by the sources
         """
 
-        bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
+        bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]], dtype = complex)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
             bPrimary[:,i] = bPrimary[:,i] + bp
@@ -277,40 +294,9 @@ class Fields_e(Fields):
             b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
         return b
 
-    def _bSecondaryDeriv_u(self, src, du_dm_v, adjoint = False):
+    def _bDeriv_u(self, src, du_dm_v, adjoint = False):
         """
-        Derivative of the secondary magnetic flux density with respect to the thing we solved for
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray du_dm_v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector
-        """
-
-        C = self._edgeCurl
-        if adjoint:
-            return - 1./(1j*omega(src.freq)) * (C.T * du_dm_v)
-        return - 1./(1j*omega(src.freq)) * (C * du_dm_v)
-
-    def _bSecondaryDeriv_m(self, src, v, adjoint = False):
-        """
-        Derivative of the secondary magnetic flux density with respect to the inversion model. 
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector
-        """
-
-        S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
-        return 1./(1j * omega(src.freq)) * S_mDeriv
-
-
-    def _bDeriv_u(self, src, du_dm_v, adjoint=False):
-        """
-        Partial derivative of the total magnetic flux density with respect to the thing we solved for
+        Derivative of the magnetic flux density with respect to the thing we solved for
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
@@ -319,21 +305,126 @@ class Fields_e(Fields):
         :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
         """
 
-        return self._bSecondaryDeriv_u(src, du_dm_v, adjoint)
+        C = self._edgeCurl
+        if adjoint:
+            return - 1./(1j*omega(src.freq)) * (C.T * du_dm_v)
+        return - 1./(1j*omega(src.freq)) * (C * du_dm_v)
 
-    def _bDeriv_m(self, src, v, adjoint=False):
+
+    def _bDeriv_m(self, src, v, adjoint = False):
         """
-        Partial derivative of the total magnetic flux density with respect to the inversion model. 
+        Derivative of the magnetic flux density with respect to the inversion model. 
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
-        :rtype: SimPEG.Utils.Zero
+        :rtype: numpy.ndarray
         :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
         """
 
-        # Assuming the primary does not depend on the model
-        return self._bSecondaryDeriv_m(src, v, adjoint)
+        S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
+        return 1./(1j * omega(src.freq)) * S_mDeriv
+
+    def _j(self,  eSolution, srcList):
+        """
+        Current density from eSolution
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: current density
+        """
+        aveE2CCV = self._aveE2CCV
+        n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        return VI * (aveE2CCV * (self._MeSigma * self._e(eSolution,srcList) ) )
+
+    def _jDeriv_u(self, src, du_dm_v, adjoint = False):
+        """
+        Derivative of the current density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the current density with respect to the field we solved for with a vector
+        """
+        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        if adjoint: 
+            return  self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
+        return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) )
+        
+
+    def _jDeriv_m(self, src, v, adjoint = False):
+        """
+        Derivative of the current density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the current density derivative with respect to the inversion model with a vector
+        """
+        e = self[src, 'e']
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        if adjoint:
+            return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
+        return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v))
+        
+
+
+    def _h(self, eSolution, srcList):
+        """
+        Magnetic field from eSolution
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: magnetic field
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
+
+    def _hDeriv_u(self, src, du_dm_v, adjoint = False):
+        """
+        Derivative of the magnetic field with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            v =  self._MfMui.T * (self._aveF2CCV.T  * (VI.T * du_dm_v))
+            return self._bDeriv_u(src, v, adjoint=adjoint)
+        return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_u(src, du_dm_v, adjoint = adjoint)))
+
+    def _hDeriv_m(self, src, v, adjoint = False):
+        """
+        Derivative of the magnetic field with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the magnetic field derivative with respect to the inversion model with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            v =  self._MfMui.T * (self._aveF2CCV.T  * (VI.T * v))
+            return self._bDeriv_m(src, v, adjoint=adjoint)
+        return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_m(src, v, adjoint = adjoint)))
+
 
 
 class Fields_b(Fields):
@@ -352,6 +443,8 @@ class Fields_b(Fields):
                     'e' : ['bSolution','E','_e'],
                     'ePrimary' : ['bSolution','E','_ePrimary'],
                     'eSecondary' : ['bSolution','E','_eSecondary'],
+                    'j' : ['bSolution','CCV','_j'],
+                    'h' : ['bSolution','CCV','_h'],
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -360,10 +453,29 @@ class Fields_b(Fields):
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
+        self._MeSigma = self.survey.prob.MeSigma
         self._MeSigmaI = self.survey.prob.MeSigmaI
         self._MfMui = self.survey.prob.MfMui
+        self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
         self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
         self._Me = self.survey.prob.Me
+        self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
+        self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
+        self._sigma = self.survey.prob.curModel.sigma
+        self._mui = self.survey.prob.curModel.mui
+        self._nC = self.survey.prob.mesh.nC
+
+
+
+    def _GLoc(self,fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
 
     def _bPrimary(self, bSolution, srcList):
         """
@@ -375,7 +487,7 @@ class Fields_b(Fields):
         :return: primary electric field as defined by the sources
         """
 
-        bPrimary = np.zeros_like(bSolution)
+        bPrimary = np.zeros([self.prob.mesh.nF,len(srcList)], dtype = complex)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
             bPrimary[:,i] = bPrimary[:,i] + bp
@@ -431,7 +543,7 @@ class Fields_b(Fields):
         :return: primary electric field as defined by the sources
         """
 
-        ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
+        ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]], dtype = complex)
         for i,src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
             ePrimary[:,i] = ePrimary[:,i] + ep
@@ -447,86 +559,139 @@ class Fields_b(Fields):
         :return: secondary electric field
         """
 
-        e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
+        e = ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
-            e[:,i] = e[:,i]+ -self._MeSigmaI * S_e
-        return e
+            e[:,i] = e[:,i] + - S_e
 
-    def _eSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
-        """
-        Derivative of the secondary electric field with respect to the thing we solved for
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary electric field with respect to the field we solved for with a vector
-        """
-
-        if not adjoint:
-            return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) )
-        else:
-            return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v))
-
-    def _eSecondaryDeriv_m(self, src, v, adjoint=False):
-        """
-        Derivative of the secondary electric field with respect to the inversion model 
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary electric field with respect to the model with a vector
-        """
-
-        bSolution = self[[src],'bSolution']
-        _,S_e = src.eval(self.prob)
-        Me = self._Me
-
-        if adjoint:
-            Me = Me.T
-
-        w = self._edgeCurl.T * (self._MfMui * bSolution)
-        w = w - Utils.mkvc(Me * S_e,2)
-
-        if not adjoint:
-            de_dm = self._MeSigmaIDeriv(w) * v
-        elif adjoint:
-            de_dm = self._MeSigmaIDeriv(w).T * v
-
-        _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
-
-        de_dm = de_dm - self._MeSigmaI * S_eDeriv
-
-        return de_dm
+        return self._MeSigmaI * e
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total electric field with respect to the thing we solved for 
+        Derivative of the electric field with respect to the thing we solved for
         
         :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray du_dm_v: vector to take product with
+        :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
         :return: product of the derivative of the electric field with respect to the field we solved for with a vector
         """
 
-        return self._eSecondaryDeriv_u(src, du_dm_v, adjoint)
+        if not adjoint:
+            return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) )
+        return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v))
+
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total electric field density with respect to the inversion model. 
+        Derivative of the electric field with respect to the inversion model 
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: product of the electric field derivative with respect to the inversion model with a vector
+        :return: product of the derivative of the electric field with respect to the model with a vector
         """
 
-        # assuming primary doesn't depend on model
-        return self._eSecondaryDeriv_m(src, v, adjoint)
+        bSolution = Utils.mkvc(self[src, 'bSolution'])
+        _,S_e = src.eval(self.prob)
+
+        w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution)
+        _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+
+
+        if adjoint:
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv
+
+    def _j(self, bSolution, srcList):
+        """
+        Secondary current density from bSolution
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary current density
+        """
+
+        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
+        
+
+    def _jDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Partial derivative of the current density with respect to the thing we 
+        solved for.
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the current density with respect to the field we solved for with a vector
+        """
+        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
+        return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) )
+
+
+    def _jDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the current density with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the current density with respect to the model with a vector
+        """
+        return Zero()
+
+    def _h(self, bSolution, srcList):
+        """
+        Magnetic field from bSolution
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: magnetic field
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
+
+    def _hDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Partial derivative of the magnetic field with respect to the thing we 
+        solved for.
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        if adjoint:
+            return self._MfMui.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v) )
+        return VI * (self._aveF2CCV * (self._MfMui * du_dm_v))
+
+    def _hDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the magnetic field with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic field with respect to the model with a vector
+        """
+        return Zero()
 
 
 class Fields_j(Fields):
@@ -545,6 +710,8 @@ class Fields_j(Fields):
                     'h' : ['jSolution','E','_h'],
                     'hPrimary' : ['jSolution','E','_hPrimary'],
                     'hSecondary' : ['jSolution','E','_hSecondary'],
+                    'e' : ['jSolution','CCV','_e'],
+                    'b' : ['jSolution','CCV','_b'],
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -553,10 +720,25 @@ class Fields_j(Fields):
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
+        self._MeMu = self.survey.prob.MeMu
         self._MeMuI = self.survey.prob.MeMuI
         self._MfRho = self.survey.prob.MfRho
         self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
-        self._Me = self.survey.prob.Me
+        self._rho = self.survey.prob.curModel.rho
+        self._mu = self.survey.prob.curModel.mui
+        self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
+        self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
+        self._nC = self.survey.prob.mesh.nC
+
+    def _GLoc(self,fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
 
     def _jPrimary(self, jSolution, srcList):
         """
@@ -568,7 +750,7 @@ class Fields_j(Fields):
         :return: primary current density as defined by the sources
         """
 
-        jPrimary = np.zeros_like(jSolution,dtype = complex)
+        jPrimary = np.zeros_like(jSolution, dtype = complex)
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
             jPrimary[:,i] = jPrimary[:,i] + jp
@@ -653,67 +835,17 @@ class Fields_j(Fields):
         :return: secondary magnetic field
         """
 
-        h =  self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) )
+        h = (self._edgeCurl.T * (self._MfRho * jSolution) )
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
             S_m,_ = src.eval(self.prob)
-            h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
-        return h
-
-
-    def _hSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
-        """
-        Derivative of the secondary magnetic field with respect to the thing we solved for
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray du_dm_v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector
-        """
-
-        if not adjoint:
-            return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
-        elif adjoint:
-            return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
-
-    def _hSecondaryDeriv_m(self, src, v, adjoint=False):
-        """
-        Derivative of the secondary magnetic field with respect to the inversion model 
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary magnetic field with respect to the model with a vector
-        """
-
-        jSolution = self[[src],'jSolution']
-        MeMuI = self._MeMuI
-        C = self._edgeCurl
-        MfRho = self._MfRho
-        MfRhoDeriv = self._MfRhoDeriv
-        Me = self._Me
-
-        if not adjoint:
-            hDeriv_m =  -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
-        elif adjoint:
-            hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
-
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
-
-        if not adjoint:
-            S_mDeriv = S_mDeriv(v)
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
-        elif adjoint:
-            S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v))
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
-        return hDeriv_m
+            h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m)
+        return self._MeMuI * h
 
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total magnetic field with respect to the thing we solved for 
+        Derivative of the magnetic field with respect to the thing we solved for
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
@@ -722,21 +854,139 @@ class Fields_j(Fields):
         :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
         """
 
-        return self._hSecondaryDeriv_u(src, du_dm_v, adjoint)
+        if adjoint:
+            return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
+        return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
+        
+
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total magnetic field density with respect to the inversion model. 
+        Derivative of the magnetic field with respect to the inversion model 
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: product of the magnetic field derivative with respect to the inversion model with a vector
+        :return: product of the derivative of the magnetic field with respect to the model with a vector
         """
 
-        # assuming the primary doesn't depend on the model
-        return self._hSecondaryDeriv_m(src, v, adjoint)
+        jSolution = Utils.mkvc(self[[src],'jSolution'])
+        MeMuI = self._MeMuI
+        C = self._edgeCurl
+        MfRho = self._MfRho
+        MfRhoDeriv = self._MfRhoDeriv
+        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+
+        if not adjoint:
+            hDeriv_m =  -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
+            S_mDeriv = S_mDeriv(v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv)
+
+        elif adjoint:
+            hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
+
+            S_mDeriv = S_mDeriv(MeMuI.T * v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
+        
+        return hDeriv_m
+
+    def _e(self, jSolution, srcList):
+        """
+        Electric field from jSolution
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: electric field
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList))) 
+
+    def _eDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Derivative of the electric field with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the electric field with respect to the field we solved for with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
+        return VI * (self._aveF2CCV * (self._MfRho *  du_dm_v))
+
+    def _eDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the electric field with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the electric field with respect to the model with a vector
+        """
+        jSolution = Utils.mkvc(self[src,'jSolution'])
+        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
+        return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) *  v))
+
+    def _b(self, jSolution, srcList):
+        """
+        Secondary magnetic flux density from jSolution
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary magnetic flux density 
+        """
+        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) 
+
+    def _bDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Derivative of the magnetic flux density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        if adjoint:
+            return -1./(1j*omega(src.freq)) * self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
+        return -1./(1j*omega(src.freq)) * VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v)))
+
+    def _bDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the magnetic flux density with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic flux density with respect to the model with a vector
+        """
+        jSolution = self[src,'jSolution']
+        n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+
+        if adjoint:
+            v = self._aveE2CCV.T * ( VI.T * v)
+            return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
+        return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
 class Fields_h(Fields):
@@ -754,7 +1004,9 @@ class Fields_h(Fields):
                     'hSecondary' : ['hSolution','E','_hSecondary'],
                     'j' : ['hSolution','F','_j'],
                     'jPrimary' : ['hSolution','F','_jPrimary'],
-                    'jSecondary' : ['hSolution','F','_jSecondary']
+                    'jSecondary' : ['hSolution','F','_jSecondary'],
+                    'e' : ['hSolution','CCV','_e'],
+                    'b' : ['hSolution','CCV','_b'],
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -763,8 +1015,25 @@ class Fields_h(Fields):
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
+        self._MeMu = self.survey.prob.MeMu
         self._MeMuI = self.survey.prob.MeMuI
         self._MfRho = self.survey.prob.MfRho
+        self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
+        self._rho = self.survey.prob.curModel.rho
+        self._mu = self.survey.prob.curModel.mui
+        self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
+        self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
+        self._nC = self.survey.prob.mesh.nC
+
+    def _GLoc(self,fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
 
     def _hPrimary(self, hSolution, srcList):
         """
@@ -841,7 +1110,7 @@ class Fields_h(Fields):
 
     def _jSecondary(self, hSolution, srcList):
         """
-        Secondary current density from eSolution
+        Secondary current density from hSolution
 
         :param numpy.ndarray hSolution: field we solved for
         :param list srcList: list of sources
@@ -855,40 +1124,9 @@ class Fields_h(Fields):
             j[:,i] = j[:,i]+ -S_e
         return j
 
-    def _jSecondaryDeriv_u(self, src, du_dm_v, adjoint=False):
-        """
-        Derivative of the secondary current density with respect to the thing we solved for
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray du_dm_v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector
-        """
-
-        if not adjoint:
-            return self._edgeCurl*du_dm_v
-        elif adjoint:
-            return self._edgeCurl.T*du_dm_v
-
-    def _jSecondaryDeriv_m(self, src, v, adjoint=False):
-        """
-        Derivative of the secondary current density with respect to the inversion model. 
-        
-        :param SimPEG.EM.FDEM.Src src: source
-        :param numpy.ndarray v: vector to take product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: product of the secondary current density derivative with respect to the inversion model with a vector
-        """
-
-        _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
-        return -S_eDeriv
-
-
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the thing we solved for
+        Derivative of the current density with respect to the thing we solved for
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
@@ -897,18 +1135,113 @@ class Fields_h(Fields):
         :return: product of the derivative of the current density with respect to the field we solved for with a vector
         """
 
-        return self._jSecondaryDeriv_u(src,du_dm_v,adjoint)
+        if not adjoint:
+            return self._edgeCurl*du_dm_v
+        elif adjoint:
+            return self._edgeCurl.T*du_dm_v
+
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the inversion model. 
+        Derivative of the current density with respect to the inversion model. 
         
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
-        :rtype: SimPEG.Utils.Zero
-        :return: product of the current density with respect to the inversion model with a vector
+        :rtype: numpy.ndarray
+        :return: product of the current density derivative with respect to the inversion model with a vector
         """
 
-        # assuming the primary does not depend on the model
-        return self._jSecondaryDeriv_m(src,v,adjoint)
+        _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+        return -S_eDeriv
+    
+    def _e(self, hSolution, srcList):
+        """
+        Electric field from hSolution
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: electric field
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList))) 
+
+    def _eDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Derivative of the electric field with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the electric field with respect to the field we solved for with a vector
+        """
+        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
+        return VI * (self._aveF2CCV * (self._MfRho *  self._edgeCurl * du_dm_v )) 
+
+    def _eDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the electric field with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the electric field derivative with respect to the inversion model with a vector
+        """
+        hSolution = Utils.mkvc(self[src,'hSolution'])
+        n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint: 
+            return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
+        return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
+
+    def _b(self, hSolution, srcList):
+        """
+        Magnetic flux density from hSolution
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: magnetic flux density 
+        """
+        h = self._h(hSolution, srcList)
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+
+        return VI * (self._aveE2CCV * (self._MeMu * h))
+
+    def _bDeriv_u(self, src, du_dm_v, adjoint=False):
+        """
+        Derivative of the magnetic flux density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray du_dm_v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
+        """
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
+        if adjoint:
+            return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
+        return VI * (self._aveE2CCV * (self._MeMu * du_dm_v))
+
+    def _bDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the magnetic flux density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
+        """
+        return Zero()
+
+
diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py
index 6193b2f4..41614b0e 100644
--- a/SimPEG/EM/FDEM/SrcFDEM.py
+++ b/SimPEG/EM/FDEM/SrcFDEM.py
@@ -335,7 +335,7 @@ class MagDipole(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
-        return h_from_b(prob,b)
+        return 1./self.mu * b 
 
     def S_m(self, prob):
         """
@@ -347,6 +347,8 @@ class MagDipole(BaseSrc):
         """
 
         b_p = self.bPrimary(prob)
+        if prob._eqLocs is 'EF':
+            b_p = prob.Me * b_p 
         return -1j*omega(self.freq)*b_p
 
     def S_e(self, prob):
@@ -449,7 +451,7 @@ class MagDipole_Bfield(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
-        return h_from_b(prob, b)
+        return 1/self.mu * b
 
     def S_m(self, prob):
         """
@@ -460,6 +462,8 @@ class MagDipole_Bfield(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
+        if prob._eqLocs is 'EF':
+            b = prob.Me * b
         return -1j*omega(self.freq)*b
 
     def S_e(self, prob):
@@ -570,6 +574,8 @@ class CircularLoop(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
+        if prob._eqLocs is 'EF':
+            b =  prob.Me *  b
         return -1j*omega(self.freq)*b
 
     def S_e(self, prob):
@@ -589,6 +595,8 @@ class CircularLoop(BaseSrc):
                 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)
@@ -596,4 +604,6 @@ class CircularLoop(BaseSrc):
 
             return -C.T * (MMui_s * self.bPrimary(prob))
 
+            
+
 
diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py
index 4d220259..d7123f71 100644
--- a/SimPEG/EM/FDEM/SurveyFDEM.py
+++ b/SimPEG/EM/FDEM/SurveyFDEM.py
@@ -3,6 +3,7 @@ from SimPEG.EM.Utils import *
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
 import SrcFDEM as Src
+from SimPEG import sp
 
 
 ####################################################
@@ -18,33 +19,33 @@ class Rx(SimPEG.Survey.BaseRx):
     """
 
     knownRxTypes = {
-                    'exr':['e', 'Ex', 'real'],
-                    'eyr':['e', 'Ey', 'real'],
-                    'ezr':['e', 'Ez', 'real'],
-                    'exi':['e', 'Ex', 'imag'],
-                    'eyi':['e', 'Ey', 'imag'],
-                    'ezi':['e', 'Ez', 'imag'],
+                    'exr':['e', 'x', 'real'],
+                    'eyr':['e', 'y', 'real'],
+                    'ezr':['e', 'z', 'real'],
+                    'exi':['e', 'x', 'imag'],
+                    'eyi':['e', 'y', 'imag'],
+                    'ezi':['e', 'z', 'imag'],
 
-                    'bxr':['b', 'Fx', 'real'],
-                    'byr':['b', 'Fy', 'real'],
-                    'bzr':['b', 'Fz', 'real'],
-                    'bxi':['b', 'Fx', 'imag'],
-                    'byi':['b', 'Fy', 'imag'],
-                    'bzi':['b', 'Fz', 'imag'],
+                    'bxr':['b', 'x', 'real'],
+                    'byr':['b', 'y', 'real'],
+                    'bzr':['b', 'z', 'real'],
+                    'bxi':['b', 'x', 'imag'],
+                    'byi':['b', 'y', 'imag'],
+                    'bzi':['b', 'z', 'imag'],
 
-                    'jxr':['j', 'Fx', 'real'],
-                    'jyr':['j', 'Fy', 'real'],
-                    'jzr':['j', 'Fz', 'real'],
-                    'jxi':['j', 'Fx', 'imag'],
-                    'jyi':['j', 'Fy', 'imag'],
-                    'jzi':['j', 'Fz', 'imag'],
+                    'jxr':['j', 'x', 'real'],
+                    'jyr':['j', 'y', 'real'],
+                    'jzr':['j', 'z', 'real'],
+                    'jxi':['j', 'x', 'imag'],
+                    'jyi':['j', 'y', 'imag'],
+                    'jzi':['j', 'z', 'imag'],
 
-                    'hxr':['h', 'Ex', 'real'],
-                    'hyr':['h', 'Ey', 'real'],
-                    'hzr':['h', 'Ez', 'real'],
-                    'hxi':['h', 'Ex', 'imag'],
-                    'hyi':['h', 'Ey', 'imag'],
-                    'hzi':['h', 'Ez', 'imag'],
+                    'hxr':['h', 'x', 'real'],
+                    'hyr':['h', 'y', 'real'],
+                    'hzr':['h', 'z', 'real'],
+                    'hxi':['h', 'x', 'imag'],
+                    'hyi':['h', 'y', 'imag'],
+                    'hzi':['h', 'z', 'imag'],
                    }
     radius = None
 
@@ -56,17 +57,16 @@ class Rx(SimPEG.Survey.BaseRx):
         """Field Type projection (e.g. e b ...)"""
         return self.knownRxTypes[self.rxType][0]
 
-    @property
-    def projGLoc(self):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        return self.knownRxTypes[self.rxType][1]
-
     @property
     def projComp(self):
         """Component projection (real/imag)"""
         return self.knownRxTypes[self.rxType][2]
 
-    def eval(self, src, mesh, f):
+    def projGLoc(self, u):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+
+    def eval(self, src, mesh, u):
         """
         Project fields to recievers to get data.
 
@@ -76,24 +76,35 @@ class Rx(SimPEG.Survey.BaseRx):
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
-        P = self.getP(mesh) # get interpolation to recievers 
-        u_part_complex = f[src, self.projField]
-        real_or_imag = self.projComp # get the real or imag component
+        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
+        # projGLoc += self.knownRxTypes[self.rxType][1]
+
+        P = self.getP(mesh, self.projGLoc(u))
+        u_part_complex = u[src, self.projField]
+        # get the real or imag component
+        real_or_imag = self.projComp
         u_part = getattr(u_part_complex, real_or_imag)
+        
         return P*u_part
 
-    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+    def evalDeriv(self, src, mesh, u, v, adjoint=False):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
         :param Source src: FDEM source
         :param Mesh mesh: mesh used
-        :param Fields f: fields object
+        :param Fields u: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
-        P = self.getP(mesh)
+
+        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
+
+        # print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc
+        # projGLoc += self.knownRxTypes[self.rxType][1]
+
+        P = self.getP(mesh, self.projGLoc(u))
 
         if not adjoint:
             Pv_complex = P * v
@@ -185,3 +196,4 @@ class Survey(SimPEG.Survey.BaseSurvey):
 
     def evalDeriv(self, u):
         raise Exception('Use Receivers to project fields deriv.')
+
diff --git a/SimPEG/EM/Utils/EMUtils.py b/SimPEG/EM/Utils/EMUtils.py
index 4a342acb..e7dbf441 100644
--- a/SimPEG/EM/Utils/EMUtils.py
+++ b/SimPEG/EM/Utils/EMUtils.py
@@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0):
     beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) 
     return alp - 1j*beta
 
-# Constitutive relations
-def e_from_j(prob,j):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MSigmaI = prob.MeSigmaI
-    elif eqLocs is 'EF':
-        MSigmaI = prob.MfRho
-    return MSigmaI*j
-
-def j_from_e(prob,e):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MSigma = prob.MeSigma
-    elif eqLocs is 'EF':
-        MSigma = prob.MfRhoI
-    return MSigma*e
-
-def b_from_h(prob,h):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MMu = prob.MfMuiI
-    elif eqLocs is 'EF':
-        MMu = prob.MeMu
-    return MMu*h
-
-def h_from_b(prob,b):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MMuI = prob.MfMui
-    elif eqLocs is 'EF':
-        MMuI = prob.MeMuI
-    return MMuI*b 
-
 
diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py
index 18dddde9..ef779042 100644
--- a/SimPEG/EM/Utils/__init__.py
+++ b/SimPEG/EM/Utils/__init__.py
@@ -1,5 +1,2 @@
-# import Sources
-# import Ana
-# import Solver
-from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
+from EMUtils import omega, k
 from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
\ No newline at end of file
diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py
index 8c703083..569f8e6d 100644
--- a/SimPEG/EM/Utils/testingUtils.py
+++ b/SimPEG/EM/Utils/testingUtils.py
@@ -4,19 +4,28 @@ 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
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = 5e-1
+
+
+def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
+    cs = 10.
+    ncx, ncy, ncz = 0, 0, 0
+    npad = 8
     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)
+    if useMu is True:
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+    else:
+        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.])
+    x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
+    XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
     Rx0 = EM.FDEM.Rx(XYZ, comp)
 
     Src = []
@@ -32,15 +41,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
             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.
+                S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
+                S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
                 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.
+                S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
+                S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
 
     if verbose:
@@ -70,6 +79,48 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
     except ImportError, e:
-        pass
+        prb.Solver = SolverLU
 
-    return prb
\ No newline at end of file
+    return prb
+
+def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False):
+
+    l2norm = lambda r: np.sqrt(r.dot(r))
+
+    prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
+    mesh = prb1.mesh
+    print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
+    
+    logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
+    mu = np.ones(mesh.nC)*MU
+
+    if addrandoms is True:
+        logsig  += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
+        mu += np.random.randn(mesh.nC)*MU*1e-1
+
+    if useMu is True:
+        m = np.r_[logsig, mu]
+    else:
+        m = logsig
+
+    survey1 = prb1.survey
+    d1 = survey1.dpred(m)
+
+    if verbose:
+        print '  Problem 1 solved'
+
+
+    prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
+
+    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(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
+    print l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol
+    return l2r < tol
diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
index 508f015c..1650b5cb 100644
--- a/SimPEG/Mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -234,6 +234,9 @@ class BaseTensorMesh(BaseMesh):
             'Fz'    -> z-component of field defined on faces
             'N'     -> scalar field defined on nodes
             'CC'    -> scalar field defined on cell centers
+            'CCVx'  -> x-component of vector field defined on cell centers
+            'CCVy'  -> y-component of vector field defined on cell centers
+            'CCVz'  -> z-component of vector field defined on cell centers
         """
         if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
             raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
@@ -257,6 +260,16 @@ class BaseTensorMesh(BaseMesh):
             Q = sp.hstack(components)
         elif locType in ['CC', 'N']:
             Q = Utils.interpmat(loc, *self.getTensor(locType))
+        elif locType in ['CCVx', 'CCVy', 'CCVz']:
+            Q = Utils.interpmat(loc, *self.getTensor('CC'))
+            Z = Utils.spzeros(loc.shape[0],self.nC)
+            if locType == 'CCVx':
+                Q = sp.hstack([Q,Z,Z])
+            elif locType == 'CCVy':
+                Q = sp.hstack([Z,Q,Z])
+            elif locType == 'CCVz':
+                Q = sp.hstack([Z,Z,Q])
+
         else:
             raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
 
diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py
index 1484b6c8..1732b2a7 100644
--- a/SimPEG/Survey.py
+++ b/SimPEG/Survey.py
@@ -35,7 +35,7 @@ class BaseRx(object):
         """Number of data in the receiver."""
         return self.locs.shape[0]
 
-    def getP(self, mesh):
+    def getP(self, mesh, projGLoc=None):
         """
             Returns the projection matrices as a
             list for all components collected by
@@ -48,7 +48,10 @@ class BaseRx(object):
         if mesh in self._Ps:
             return self._Ps[mesh]
 
-        P = mesh.getInterpolationMat(self.locs, self.projGLoc)
+        if projGLoc is None:
+            projGLoc = self.projGLoc
+
+        P = mesh.getInterpolationMat(self.locs, projGLoc)
         if self.storeProjections:
             self._Ps[mesh] = P
         return P
diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py
index 437f3708..da446675 100644
--- a/tests/em/fdem/forward/test_FDEM_forward.py
+++ b/tests/em/fdem/forward/test_FDEM_forward.py
@@ -3,125 +3,75 @@ from SimPEG import *
 from SimPEG import EM
 import sys
 from scipy.constants import mu_0
-from SimPEG.EM.Utils.testingUtils import getFDEMProblem
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
 
 testEB = True
 testHJ = True
-
+testEJ = True
+testBH = 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
+TOLEBHJ = 1e-5
+TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
+#TODO: choose better testing parameters to lower this 
 
 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'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose))
         def test_EB_CrossCheck_eyr_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'eyr'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose))
         def test_EB_CrossCheck_ezr_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'ezr'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose))
         def test_EB_CrossCheck_exi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'exi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose))
         def test_EB_CrossCheck_eyi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'eyi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose))
         def test_EB_CrossCheck_ezi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'ezi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose))
 
         def test_EB_CrossCheck_bxr_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'bxr'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose))
         def test_EB_CrossCheck_byr_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'byr'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose))
         def test_EB_CrossCheck_bzr_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'bzr'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose))
         def test_EB_CrossCheck_bxi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'bxi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose))
         def test_EB_CrossCheck_byi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'byi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose))
         def test_EB_CrossCheck_bzi_Eform(self):
-            self.assertTrue(crossCheckTest('e', 'bzi'))
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose))
 
     if testHJ:
         def test_HJ_CrossCheck_jxr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jxr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxr', verbose=verbose))
         def test_HJ_CrossCheck_jyr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jyr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyr', verbose=verbose))
         def test_HJ_CrossCheck_jzr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jzr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzr', verbose=verbose))
         def test_HJ_CrossCheck_jxi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jxi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxi', verbose=verbose))
         def test_HJ_CrossCheck_jyi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jyi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyi', verbose=verbose))
         def test_HJ_CrossCheck_jzi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'jzi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzi', verbose=verbose))
 
         def test_HJ_CrossCheck_hxr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hxr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxr', verbose=verbose))
         def test_HJ_CrossCheck_hyr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hyr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyr', verbose=verbose))
         def test_HJ_CrossCheck_hzr_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hzr'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzr', verbose=verbose))
         def test_HJ_CrossCheck_hxi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hxi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxi', verbose=verbose))
         def test_HJ_CrossCheck_hyi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hyi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyi', verbose=verbose))
         def test_HJ_CrossCheck_hzi_Jform(self):
-            self.assertTrue(crossCheckTest('j', 'hzi'))
+            self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzi', verbose=verbose))
 
 if __name__ == '__main__':
     unittest.main()
\ No newline at end of file
diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py
new file mode 100644
index 00000000..e6319dfc
--- /dev/null
+++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py
@@ -0,0 +1,125 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
+
+testEJ = True
+testBH = True
+
+TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
+#TODO: choose better testing parameters to lower this 
+
+SrcList = ['RawVec', 'MagDipole', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
+
+
+class FDEM_CrossCheck(unittest.TestCase):
+    if testEJ:
+        def test_EJ_CrossCheck_jxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', TOL=TOLEJHB))
+        def test_EJ_CrossCheck_jyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_jzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_jxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_jyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_jzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi',  TOL=TOLEJHB))
+
+        def test_EJ_CrossCheck_exr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_eyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_ezr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_exi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_eyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_ezi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi',  TOL=TOLEJHB))
+
+        def test_EJ_CrossCheck_bxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_byr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_bzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_bxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_byi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_bzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi',  TOL=TOLEJHB))
+
+        def test_EJ_CrossCheck_hxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_hyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_hzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_hxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_hyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi',  TOL=TOLEJHB))
+        def test_EJ_CrossCheck_hzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi',  TOL=TOLEJHB))
+
+    if testBH:
+        def test_HB_CrossCheck_jxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxr', TOL=TOLEJHB))
+        def test_HB_CrossCheck_jyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_jzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_jxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_jyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_jzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzi',  TOL=TOLEJHB))
+
+        def test_HB_CrossCheck_exr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_eyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_ezr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_exi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_eyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_ezi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezi',  TOL=TOLEJHB))
+
+        def test_HB_CrossCheck_bxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_byr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_bzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_bxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_byi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_bzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzi',  TOL=TOLEJHB))
+
+        def test_HB_CrossCheck_hxr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_hyr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_hzr_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzr',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_hxi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_hyi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyi',  TOL=TOLEJHB))
+        def test_HB_CrossCheck_hzi_Jform(self):
+            self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzi',  TOL=TOLEJHB))
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py
new file mode 100644
index 00000000..545a5014
--- /dev/null
+++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py
@@ -0,0 +1,128 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest
+
+testEB = True
+testHJ = True
+testEJ = True
+testBH = True
+verbose = False
+
+TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
+#TODO: choose better testing parameters to lower this 
+
+SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
+
+
+class FDEM_CrossCheck(unittest.TestCase):
+    if testBH:
+        def test_BH_CrossCheck_jxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_exr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_eyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_ezr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_exi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_eyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_ezi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_bxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_byr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_byi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_hxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
+
+    if testBH:
+        def test_BH_CrossCheck_jxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_jzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_exr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_eyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_ezr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_exi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_eyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_ezi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_bxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_byr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_byi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_bzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB))
+
+        def test_BH_CrossCheck_hxr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hyr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hzr(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hxi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hyi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB))
+        def test_BH_CrossCheck_hzi(self):
+            self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py
similarity index 57%
rename from tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py
rename to tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py
index f77f2131..25762368 100644
--- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py
+++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py
@@ -5,8 +5,8 @@ import sys
 from scipy.constants import mu_0
 from SimPEG.EM.Utils.testingUtils import getFDEMProblem
 
-testEB = True
-testHJ = True
+testE = True
+testB = True
 
 verbose = False
 
@@ -17,10 +17,10 @@ MU = mu_0
 freq = 1e-1
 addrandoms = True
 
-SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
 
 def adjointTest(fdemType, comp):
-    prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
+    prb = getFDEMProblem(fdemType, comp, SrcList, freq)
     print 'Adjoint %s formulation - %s' % (fdemType, comp)
 
     m  = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
@@ -45,7 +45,7 @@ def adjointTest(fdemType, comp):
     return np.abs(vJw - wJtv) < tol
 
 class FDEM_AdjointTests(unittest.TestCase):
-    if testEB:
+    if testE:
         def test_Jtvec_adjointTest_exr_Eform(self):
             self.assertTrue(adjointTest('e', 'exr'))
         def test_Jtvec_adjointTest_eyr_Eform(self):
@@ -72,6 +72,33 @@ class FDEM_AdjointTests(unittest.TestCase):
         def test_Jtvec_adjointTest_bzi_Eform(self):
             self.assertTrue(adjointTest('e', 'bzi'))
 
+        def test_Jtvec_adjointTest_jxr_Eform(self):
+            self.assertTrue(adjointTest('e', 'jxr'))
+        def test_Jtvec_adjointTest_jyr_Eform(self):
+            self.assertTrue(adjointTest('e', 'jyr'))
+        def test_Jtvec_adjointTest_jzr_Eform(self):
+            self.assertTrue(adjointTest('e', 'jzr'))
+        def test_Jtvec_adjointTest_jxi_Eform(self):
+            self.assertTrue(adjointTest('e', 'jxi'))
+        def test_Jtvec_adjointTest_jyi_Eform(self):
+            self.assertTrue(adjointTest('e', 'jyi'))
+        def test_Jtvec_adjointTest_jzi_Eform(self):
+            self.assertTrue(adjointTest('e', 'jzi'))
+
+        def test_Jtvec_adjointTest_hxr_Eform(self):
+            self.assertTrue(adjointTest('e', 'hxr'))
+        def test_Jtvec_adjointTest_hyr_Eform(self):
+            self.assertTrue(adjointTest('e', 'hyr'))
+        def test_Jtvec_adjointTest_hzr_Eform(self):
+            self.assertTrue(adjointTest('e', 'hzr'))
+        def test_Jtvec_adjointTest_hxi_Eform(self):
+            self.assertTrue(adjointTest('e', 'hxi'))
+        def test_Jtvec_adjointTest_hyi_Eform(self):
+            self.assertTrue(adjointTest('e', 'hyi'))
+        def test_Jtvec_adjointTest_hzi_Eform(self):
+            self.assertTrue(adjointTest('e', 'hzi'))
+
+    if testB:
         def test_Jtvec_adjointTest_exr_Bform(self):
             self.assertTrue(adjointTest('b', 'exr'))
         def test_Jtvec_adjointTest_eyr_Bform(self):
@@ -84,6 +111,7 @@ class FDEM_AdjointTests(unittest.TestCase):
             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):
@@ -97,59 +125,31 @@ class FDEM_AdjointTests(unittest.TestCase):
         def test_Jtvec_adjointTest_bzi_Bform(self):
             self.assertTrue(adjointTest('b', 'bzi'))
 
+        def test_Jtvec_adjointTest_jxr_Bform(self):
+            self.assertTrue(adjointTest('b', 'jxr'))
+        def test_Jtvec_adjointTest_jyr_Bform(self):
+            self.assertTrue(adjointTest('b', 'jyr'))
+        def test_Jtvec_adjointTest_jzr_Bform(self):
+            self.assertTrue(adjointTest('b', 'jzr'))
+        def test_Jtvec_adjointTest_jxi_Bform(self):
+            self.assertTrue(adjointTest('b', 'jxi'))
+        def test_Jtvec_adjointTest_jyi_Bform(self):
+            self.assertTrue(adjointTest('b', 'jyi'))
+        def test_Jtvec_adjointTest_jzi_Bform(self):
+            self.assertTrue(adjointTest('b', 'jzi'))
 
-    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'))
+        def test_Jtvec_adjointTest_hxr_Bform(self):
+            self.assertTrue(adjointTest('b', 'hxr'))
+        def test_Jtvec_adjointTest_hyr_Bform(self):
+            self.assertTrue(adjointTest('b', 'hyr'))
+        def test_Jtvec_adjointTest_hzr_Bform(self):
+            self.assertTrue(adjointTest('b', 'hzr'))
+        def test_Jtvec_adjointTest_hxi_Bform(self):
+            self.assertTrue(adjointTest('b', 'hxi'))
+        def test_Jtvec_adjointTest_hyi_Bform(self):
+            self.assertTrue(adjointTest('b', 'hyi'))
+        def test_Jtvec_adjointTest_hzi_Bform(self):
+            self.assertTrue(adjointTest('b', 'hzi'))
 
 
 if __name__ == '__main__':
diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py
new file mode 100644
index 00000000..c3fb3d37
--- /dev/null
+++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py
@@ -0,0 +1,155 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+import sys
+from scipy.constants import mu_0
+from SimPEG.EM.Utils.testingUtils import getFDEMProblem
+
+testJ = True
+testH = 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'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+
+def adjointTest(fdemType, comp):
+    prb = getFDEMProblem(fdemType, comp, SrcList, 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
+    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 testJ:
+        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_exr_Jform(self):
+            self.assertTrue(adjointTest('j', 'exr'))
+        def test_Jtvec_adjointTest_eyr_Jform(self):
+            self.assertTrue(adjointTest('j', 'eyr'))
+        def test_Jtvec_adjointTest_ezr_Jform(self):
+            self.assertTrue(adjointTest('j', 'ezr'))
+        def test_Jtvec_adjointTest_exi_Jform(self):
+            self.assertTrue(adjointTest('j', 'exi'))
+        def test_Jtvec_adjointTest_eyi_Jform(self):
+            self.assertTrue(adjointTest('j', 'eyi'))
+        def test_Jtvec_adjointTest_ezi_Jform(self):
+            self.assertTrue(adjointTest('j', 'ezi'))
+
+        def test_Jtvec_adjointTest_bxr_Jform(self):
+            self.assertTrue(adjointTest('j', 'bxr'))
+        def test_Jtvec_adjointTest_byr_Jform(self):
+            self.assertTrue(adjointTest('j', 'byr'))
+        def test_Jtvec_adjointTest_bzr_Jform(self):
+            self.assertTrue(adjointTest('j', 'bzr'))
+        def test_Jtvec_adjointTest_bxi_Jform(self):
+            self.assertTrue(adjointTest('j', 'bxi'))
+        def test_Jtvec_adjointTest_byi_Jform(self):
+            self.assertTrue(adjointTest('j', 'byi'))
+        def test_Jtvec_adjointTest_bzi_Jform(self):
+            self.assertTrue(adjointTest('j', 'bzi'))
+
+    if testH: 
+        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_jxr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jxr'))
+        def test_Jtvec_adjointTest_jyr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jyr'))
+        def test_Jtvec_adjointTest_jzr_Hform(self):
+            self.assertTrue(adjointTest('h', 'jzr'))
+        def test_Jtvec_adjointTest_jxi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jxi'))
+        def test_Jtvec_adjointTest_jyi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jyi'))
+        def test_Jtvec_adjointTest_jzi_Hform(self):
+            self.assertTrue(adjointTest('h', 'jzi'))
+
+        def test_Jtvec_adjointTest_exr_Hform(self):
+            self.assertTrue(adjointTest('h', 'exr'))
+        def test_Jtvec_adjointTest_eyr_Hform(self):
+            self.assertTrue(adjointTest('h', 'eyr'))
+        def test_Jtvec_adjointTest_ezr_Hform(self):
+            self.assertTrue(adjointTest('h', 'ezr'))
+        def test_Jtvec_adjointTest_exi_Hform(self):
+            self.assertTrue(adjointTest('h', 'exi'))
+        def test_Jtvec_adjointTest_eyi_Hform(self):
+            self.assertTrue(adjointTest('h', 'eyi'))
+        def test_Jtvec_adjointTest_ezi_Hform(self):
+            self.assertTrue(adjointTest('h', 'ezi'))
+
+        def test_Jtvec_adjointTest_bxr_Hform(self):
+            self.assertTrue(adjointTest('h', 'bxr'))
+        def test_Jtvec_adjointTest_byr_Hform(self):
+            self.assertTrue(adjointTest('h', 'byr'))
+        def test_Jtvec_adjointTest_bzr_Hform(self):
+            self.assertTrue(adjointTest('h', 'bzr'))
+        def test_Jtvec_adjointTest_bxi_Hform(self):
+            self.assertTrue(adjointTest('h', 'bxi'))
+        def test_Jtvec_adjointTest_byi_Hform(self):
+            self.assertTrue(adjointTest('h', 'byi'))
+        def test_Jtvec_adjointTest_bzi_Hform(self):
+            self.assertTrue(adjointTest('h', 'bzi'))
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
index d3bcb218..0a2e8b82 100644
--- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
+++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
@@ -5,9 +5,11 @@ import sys
 from scipy.constants import mu_0
 from SimPEG.EM.Utils.testingUtils import getFDEMProblem
 
-testDerivs = True
-testEB = True
-testHJ = True
+
+testE = True
+testB = True
+testH = True
+testJ = True
 
 verbose = False
 
@@ -18,12 +20,12 @@ MU = mu_0
 freq = 1e-1
 addrandoms = True
 
-SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
 
 
 def derivTest(fdemType, comp):
 
-    prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
+    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)
@@ -32,9 +34,6 @@ def derivTest(fdemType, comp):
         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)
@@ -43,7 +42,7 @@ def derivTest(fdemType, comp):
 
 class FDEM_DerivTests(unittest.TestCase):
 
-    if testEB:
+    if testE:
         def test_Jvec_exr_Eform(self):
             self.assertTrue(derivTest('e', 'exr'))
         def test_Jvec_eyr_Eform(self):
@@ -70,6 +69,33 @@ class FDEM_DerivTests(unittest.TestCase):
         def test_Jvec_bzi_Eform(self):
             self.assertTrue(derivTest('e', 'bzi'))
 
+        def test_Jvec_exr_Eform(self):
+            self.assertTrue(derivTest('e', 'jxr'))
+        def test_Jvec_eyr_Eform(self):
+            self.assertTrue(derivTest('e', 'jyr'))
+        def test_Jvec_ezr_Eform(self):
+            self.assertTrue(derivTest('e', 'jzr'))
+        def test_Jvec_exi_Eform(self):
+            self.assertTrue(derivTest('e', 'jxi'))
+        def test_Jvec_eyi_Eform(self):
+            self.assertTrue(derivTest('e', 'jyi'))
+        def test_Jvec_ezi_Eform(self):
+            self.assertTrue(derivTest('e', 'jzi'))
+
+        def test_Jvec_bxr_Eform(self):
+            self.assertTrue(derivTest('e', 'hxr'))
+        def test_Jvec_byr_Eform(self):
+            self.assertTrue(derivTest('e', 'hyr'))
+        def test_Jvec_bzr_Eform(self):
+            self.assertTrue(derivTest('e', 'hzr'))
+        def test_Jvec_bxi_Eform(self):
+            self.assertTrue(derivTest('e', 'hxi'))
+        def test_Jvec_byi_Eform(self):
+            self.assertTrue(derivTest('e', 'hyi'))
+        def test_Jvec_bzi_Eform(self):
+            self.assertTrue(derivTest('e', 'hzi'))
+
+    if testB:
         def test_Jvec_exr_Bform(self):
             self.assertTrue(derivTest('b', 'exr'))
         def test_Jvec_eyr_Bform(self):
@@ -96,7 +122,33 @@ class FDEM_DerivTests(unittest.TestCase):
         def test_Jvec_bzi_Bform(self):
             self.assertTrue(derivTest('b', 'bzi'))
 
-    if testHJ:
+        def test_Jvec_jxr_Bform(self):
+            self.assertTrue(derivTest('b', 'jxr'))
+        def test_Jvec_jyr_Bform(self):
+            self.assertTrue(derivTest('b', 'jyr'))
+        def test_Jvec_jzr_Bform(self):
+            self.assertTrue(derivTest('b', 'jzr'))
+        def test_Jvec_jxi_Bform(self):
+            self.assertTrue(derivTest('b', 'jxi'))
+        def test_Jvec_jyi_Bform(self):
+            self.assertTrue(derivTest('b', 'jyi'))
+        def test_Jvec_jzi_Bform(self):
+            self.assertTrue(derivTest('b', 'jzi'))
+
+        def test_Jvec_hxr_Bform(self):
+            self.assertTrue(derivTest('b', 'hxr'))
+        def test_Jvec_hyr_Bform(self):
+            self.assertTrue(derivTest('b', 'hyr'))
+        def test_Jvec_hzr_Bform(self):
+            self.assertTrue(derivTest('b', 'hzr'))
+        def test_Jvec_hxi_Bform(self):
+            self.assertTrue(derivTest('b', 'hxi'))
+        def test_Jvec_hyi_Bform(self):
+            self.assertTrue(derivTest('b', 'hyi'))
+        def test_Jvec_hzi_Bform(self):
+            self.assertTrue(derivTest('b', 'hzi'))
+
+    if testJ:
         def test_Jvec_jxr_Jform(self):
             self.assertTrue(derivTest('j', 'jxr'))
         def test_Jvec_jyr_Jform(self):
@@ -123,6 +175,34 @@ class FDEM_DerivTests(unittest.TestCase):
         def test_Jvec_hzi_Jform(self):
             self.assertTrue(derivTest('j', 'hzi'))
 
+        def test_Jvec_exr_Jform(self):
+            self.assertTrue(derivTest('j', 'exr'))
+        def test_Jvec_eyr_Jform(self):
+            self.assertTrue(derivTest('j', 'eyr'))
+        def test_Jvec_ezr_Jform(self):
+            self.assertTrue(derivTest('j', 'ezr'))
+        def test_Jvec_exi_Jform(self):
+            self.assertTrue(derivTest('j', 'exi'))
+        def test_Jvec_eyi_Jform(self):
+            self.assertTrue(derivTest('j', 'eyi'))
+        def test_Jvec_ezi_Jform(self):
+            self.assertTrue(derivTest('j', 'ezi'))
+
+        def test_Jvec_bxr_Jform(self):
+            self.assertTrue(derivTest('j', 'bxr'))
+        def test_Jvec_byr_Jform(self):
+            self.assertTrue(derivTest('j', 'byr'))
+        def test_Jvec_bzr_Jform(self):
+            self.assertTrue(derivTest('j', 'bzr'))
+        def test_Jvec_bxi_Jform(self):
+            self.assertTrue(derivTest('j', 'bxi'))
+        def test_Jvec_byi_Jform(self):
+            self.assertTrue(derivTest('j', 'byi'))
+        def test_Jvec_bzi_Jform(self):
+            self.assertTrue(derivTest('j', 'bzi'))
+
+
+    if testH:
         def test_Jvec_hxr_Hform(self):
             self.assertTrue(derivTest('h', 'hxr'))
         def test_Jvec_hyr_Hform(self):
@@ -149,6 +229,32 @@ class FDEM_DerivTests(unittest.TestCase):
         def test_Jvec_hzi_Hform(self):
             self.assertTrue(derivTest('h', 'jzi'))
 
+        def test_Jvec_exr_Hform(self):
+            self.assertTrue(derivTest('h', 'exr'))
+        def test_Jvec_eyr_Hform(self):
+            self.assertTrue(derivTest('h', 'eyr'))
+        def test_Jvec_ezr_Hform(self):
+            self.assertTrue(derivTest('h', 'ezr'))
+        def test_Jvec_exi_Hform(self):
+            self.assertTrue(derivTest('h', 'exi'))
+        def test_Jvec_eyi_Hform(self):
+            self.assertTrue(derivTest('h', 'eyi'))
+        def test_Jvec_ezi_Hform(self):
+            self.assertTrue(derivTest('h', 'ezi'))
+
+        def test_Jvec_bxr_Hform(self):
+            self.assertTrue(derivTest('h', 'bxr'))
+        def test_Jvec_byr_Hform(self):
+            self.assertTrue(derivTest('h', 'byr'))
+        def test_Jvec_bzr_Hform(self):
+            self.assertTrue(derivTest('h', 'bzr'))
+        def test_Jvec_bxi_Hform(self):
+            self.assertTrue(derivTest('h', 'bxi'))
+        def test_Jvec_byi_Hform(self):
+            self.assertTrue(derivTest('h', 'byi'))
+        def test_Jvec_bzi_Hform(self):
+            self.assertTrue(derivTest('h', 'bzi'))
+
 
 if __name__ == '__main__':
     unittest.main()