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/Utils/EMUtils.py b/SimPEG/EM/Utils/EMUtils.py
index 4a342acb..660ef117 100644
--- a/SimPEG/EM/Utils/EMUtils.py
+++ b/SimPEG/EM/Utils/EMUtils.py
@@ -13,37 +13,37 @@ 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
+# # 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 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 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 
+# 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..a60d9701 100644
--- a/SimPEG/EM/Utils/__init__.py
+++ b/SimPEG/EM/Utils/__init__.py
@@ -1,5 +1,5 @@
 # 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 ae2b7321..569f8e6d 100644
--- a/SimPEG/EM/Utils/testingUtils.py
+++ b/SimPEG/EM/Utils/testingUtils.py
@@ -8,10 +8,9 @@ FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
 CONDUCTIVITY = 1e1
 MU = mu_0
 freq = 5e-1
-addrandoms = False
 
 
-def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
+def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     cs = 10.
     ncx, ncy, ncz = 0, 0, 0
     npad = 8
@@ -20,9 +19,12 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
     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(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) #don't sample right by the source
+    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)
 
@@ -81,22 +83,26 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
 
     return prb
 
-def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False):
+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, verbose)
+    prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
     print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
-    m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
-    mu = np.log(np.ones(mesh.nC)*MU)
+    
+    logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
+    mu = 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
+        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
 
-    # prb1.PropMap.PropModel.mu = mu
-    # prb1.PropMap.PropModel.mui = 1./mu
     survey1 = prb1.survey
     d1 = survey1.dpred(m)
 
@@ -104,9 +110,8 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False)
         print '  Problem 1 solved'
 
 
-    prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose)
+    prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
 
-    # prb2.mu = mu
     survey2 = prb2.survey
     d2 = survey2.dpred(m)
 
@@ -116,6 +121,6 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False)
     r = d2-d1
     l2r = l2norm(r)
 
-    tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
+    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/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py
index d4218b86..e6319dfc 100644
--- a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py
+++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py
@@ -11,7 +11,7 @@ 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_Bfield', 'MagDipole', 'CircularLoop']
+SrcList = ['RawVec', 'MagDipole', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
 
 
 class FDEM_CrossCheck(unittest.TestCase):