diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py
index 36a2de63..9abb0a15 100644
--- a/simpegEM/Analytics/FDEM.py
+++ b/simpegEM/Analytics/FDEM.py
@@ -5,7 +5,8 @@ from scipy.special import erf
 import matplotlib.pyplot as plt
 from SimPEG import Utils
 
-def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
+
+def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0):
     """
     4.56 in Ward and Hohmann
 
@@ -25,7 +26,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
 
     """
     r = np.abs(r)
-    k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma)
+    k = np.sqrt(-1j*2.*np.pi*freq*mu*sigma)
 
     m = 1
     front = m / (2. * np.pi * (k**2) * (r**5) )
@@ -41,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
 
     return hz
 
-def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
+def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
     """
     Analytical solution for a dipole in a whole-space.
 
@@ -75,10 +76,10 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
     dz = XYZ[:,2]-srcLoc[2]
 
     r  = np.sqrt( dx**2. + dy**2. + dz**2.)
-    k  = np.sqrt( -1j*2.*np.pi*f*mu_0*sig )
+    k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
     kr = k*r
 
-    front = m / (4.*pi * r**3.) * np.exp(-1j*kr)
+    front = moment / (4.*pi * r**3.) * np.exp(-1j*kr)
     mid   = -kr**2. + 3.*1j*kr + 3.
 
     if orientation.upper() == 'X':
@@ -96,9 +97,9 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
         Hy = front*( (dy*dz/r**2.) * mid )
         Hz = front*( (dz/r)**2. * mid + (kr**2. - 1j*kr - 1.) )
 
-    Bx = mu_0*Hx
-    By = mu_0*Hy
-    Bz = mu_0*Hz
+    Bx = mu*Hx
+    By = mu*Hy
+    Bz = mu*Hz
 
     if Bx.ndim is 1:
         Bx = Utils.mkvc(Bx,2)
@@ -112,7 +113,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
     return Bx, By, Bz
 
 
-def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
+def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0):
     XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
 
     dx = XYZ[:,0]-srcLoc[0]
@@ -120,21 +121,33 @@ def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
     dz = XYZ[:,2]-srcLoc[2]
 
     r  = np.sqrt( dx**2. + dy**2. + dz**2.)
-    k  = np.sqrt( -1j*2.*np.pi*f*mu_0*sig )
+    k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
     kr = k*r 
 
-    front = moment / (4. * np.pi * sig * r**3) * exp(-1j*k*r)
+    front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*k*r)
     mid   = -k**2 * r**2 + 3*1j*k*r + 3 
 
-    Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r))
-    Ey = front*(dx*dy  / r**2)*mid
-    Ez = front*(dx*dz  / r**2)*mid
+    # Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r))
+    # Ey = front*(dx*dy  / r**2)*mid
+    # Ez = front*(dx*dz  / r**2)*mid
 
     if orientation.upper() == 'X':
+        Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ey = front*(dx*dy  / r**2)*mid
+        Ez = front*(dx*dz  / r**2)*mid   
         return Ex, Ey, Ez
 
     elif orientation.upper() == 'Y':
-        return Ez, Ex, Ey 
+        #  x--> y, y--> z, z-->x
+        Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ez = front*(dy*dz  / r**2)*mid
+        Ex = front*(dy*dx  / r**2)*mid   
+        return Ex, Ey, Ez
 
     elif orientation.upper() == 'Z':
-        return Ey, Ez, Ex 
\ No newline at end of file
+        # x --> z, y --> x, z --> y 
+        Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ex = front*(dz*dx  / r**2)*mid
+        Ey = front*(dz*dy  / r**2)*mid   
+        return Ex, Ey, Ez
+        # return Ey, Ez, Ex 
\ No newline at end of file
diff --git a/simpegEM/Base.py b/simpegEM/Base.py
index 4205c518..5a17c69b 100644
--- a/simpegEM/Base.py
+++ b/simpegEM/Base.py
@@ -141,5 +141,5 @@ class BaseEMProblem(Problem.BaseProblem):
 
     # TODO: This isn't going to work yet
     # TODO: This should take a vector
-    def dMfRhoIDeriv(self,u):
+    def MfRhoIDeriv(self,u):
         return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index 5e465f36..dde77f15 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -19,7 +19,7 @@ class BaseFDEMProblem(BaseEMProblem):
     surveyPair = SurveyFDEM
     fieldsPair = FieldsFDEM
 
-    def fields(self, m):
+    def fields(self, m=None):
         self.curModel = m
         F = self.fieldsPair(self.mesh, self.survey)
 
@@ -151,10 +151,6 @@ class BaseFDEMProblem(BaseEMProblem):
 
         return S_m, S_e
 
-    def getSourceTermDeriv(self,freq,m,v,u=None,adjoint=False):
-        raise NotImplementedError('getSourceTermDeriv not implemented yet')
-        return None, None
-
 
 ##########################################################################################
 ################################ E-B Formulation #########################################
@@ -321,14 +317,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
 
     def getRHSDeriv_m(self, src, v, adjoint=False):
         C = self.mesh.edgeCurl
-        S_m, S_e = self.getSourceTerm(src.freq)
+        S_m, S_e = src.eval(self)
         MfMui = self.MfMui
 
         if self._makeASymmetric and adjoint:
             v = self.MfMui * v
 
         if S_e is not None:
-            MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
+            MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc(S_e))
             if not adjoint:
                 RHSderiv = C * (MeSigmaIDeriv * v)
             elif adjoint:
@@ -577,7 +573,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
         return RHS
 
     def getRHSDeriv_m(self, src, v, adjoint=False):
-        _, S_e = self.getSourceTerm(src.freq)
+        _, S_e = src.eval(self)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
         MfRhoDeriv = self.MfRhoDeriv(S_e)
diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py
index 5838ddbe..46b67fa7 100644
--- a/simpegEM/FDEM/FieldsFDEM.py
+++ b/simpegEM/FDEM/FieldsFDEM.py
@@ -7,7 +7,6 @@ class FieldsFDEM(Problem.Fields):
     knownFields = {}
     dtype = complex
 
-
 class FieldsFDEM_e(FieldsFDEM):
     knownFields = {'eSolution':'E'}
     aliasFields = {
@@ -23,6 +22,7 @@ class FieldsFDEM_e(FieldsFDEM):
         FieldsFDEM.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
+        self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
 
     # def getDeriv_u(self, fieldsList, src, v, adjoint=False):
@@ -32,7 +32,7 @@ class FieldsFDEM_e(FieldsFDEM):
     def _ePrimary(self, eSolution, srcList):
         ePrimary = np.zeros_like(eSolution)
         for i, src in enumerate(srcList):
-            ep = src.ePrimary(self.survey.prob)
+            ep = src.ePrimary(self.prob)
             if ep is not None:
                 ePrimary[:,i] = ep     
         return ePrimary
@@ -53,7 +53,7 @@ class FieldsFDEM_e(FieldsFDEM):
     def _bPrimary(self, eSolution, srcList):
         bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
-            bp = src.bPrimary(self.survey.prob)
+            bp = src.bPrimary(self.prob)
             if bp is not None:
                 bPrimary[:,i] += bp
         return bPrimary
@@ -63,7 +63,7 @@ class FieldsFDEM_e(FieldsFDEM):
         b = (C * eSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            S_m, _ = src.eval(self.survey.prob)
+            S_m, _ = src.eval(self.prob)
             if S_m is not None:
                 b[:,i] += 1./(1j*omega(src.freq)) * S_m
         return b
@@ -75,7 +75,7 @@ class FieldsFDEM_e(FieldsFDEM):
         return - 1./(1j*omega(src.freq)) * (C * v)
 
     def _bSecondaryDeriv_m(self, src, v, adjoint = False):
-        S_mDeriv, _ = src.evalDeriv(self.survey.prob, adjoint)
+        S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
         S_mDeriv = S_mDeriv(v)
         if S_mDeriv is not None:
             return 1./(1j * omega(src.freq)) * S_mDeriv
@@ -108,6 +108,7 @@ class FieldsFDEM_b(FieldsFDEM):
         FieldsFDEM.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
+        self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
         self._MeSigmaI = self.survey.prob.MeSigmaI
         self._MfMui = self.survey.prob.MfMui
@@ -116,7 +117,7 @@ class FieldsFDEM_b(FieldsFDEM):
     def _bPrimary(self, bSolution, srcList):
         bPrimary = np.zeros_like(bSolution)
         for i, src in enumerate(srcList):
-            bp = src.bPrimary(self.survey.prob)
+            bp = src.bPrimary(self.prob)
             if bp is not None:
                 bPrimary[:,i] = bp
         return bPrimary
@@ -137,7 +138,7 @@ class FieldsFDEM_b(FieldsFDEM):
     def _ePrimary(self, bSolution, srcList):
         ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
         for i,src in enumerate(srcList):
-            ep = src.ePrimary(self.survey.prob)
+            ep = src.ePrimary(self.prob)
             if ep is not None:
                 ePrimary[:,i] = ep
         return ePrimary
@@ -145,9 +146,9 @@ class FieldsFDEM_b(FieldsFDEM):
     def _eSecondary(self, bSolution, srcList):
         e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList): 
-            _,S_e = src.eval(self.survey.prob)
+            _,S_e = src.eval(self.prob)
             if S_e is not None:
-                e += -self._MeSigmaI*S_e
+                e[:,i] += -self._MeSigmaI*S_e
         return e
 
     def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -158,18 +159,18 @@ class FieldsFDEM_b(FieldsFDEM):
 
     def _eSecondaryDeriv_m(self, src, v, adjoint=False):
         bSolution = self[[src],'bSolution']
-        _,S_e = src.eval(self.survey.prob)
+        _,S_e = src.eval(self.prob)
 
         w = self._edgeCurl.T * (self._MfMui * bSolution)
         if S_e is not None:
-            w += -S_e
+            w += -Utils.mkvc(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.survey.prob, adjoint)
+        _, S_eDeriv = src.evalDeriv(self.prob, adjoint)
         Se_Deriv = S_eDeriv(v)
 
         if Se_Deriv is not None:
@@ -203,6 +204,7 @@ class FieldsFDEM_j(FieldsFDEM):
         FieldsFDEM.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
+        self.prob = self.survey.prob 
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
         self._MeMuI = self.survey.prob.MeMuI
         self._MfRho = self.survey.prob.MfRho
@@ -211,7 +213,7 @@ class FieldsFDEM_j(FieldsFDEM):
     def _jPrimary(self, jSolution, srcList):
         jPrimary = np.zeros_like(jSolution,dtype = complex)
         for i, src in enumerate(srcList):
-            jp = src.jPrimary(self.survey.prob) 
+            jp = src.jPrimary(self.prob) 
             if jp is not None:
                 jPrimary[:,i] += jp
         return jPrimary
@@ -232,7 +234,7 @@ class FieldsFDEM_j(FieldsFDEM):
     def _hPrimary(self, jSolution, srcList):
         hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
-            hp = src.hPrimary(self.survey.prob)
+            hp = src.hPrimary(self.prob)
             if hp is not None:
                 hPrimary[:,i] = hp 
         return hPrimary
@@ -244,7 +246,7 @@ class FieldsFDEM_j(FieldsFDEM):
         h =  MeMuI * (C.T * (MfRho * jSolution) ) 
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
-            S_m,_ = src.eval(self.survey.prob)
+            S_m,_ = src.eval(self.prob)
             if S_m is not None:
                 h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
         return h
@@ -270,7 +272,7 @@ class FieldsFDEM_j(FieldsFDEM):
         elif adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) 
 
-        S_mDeriv,_ = src.evalDeriv(self.survey.prob, adjoint)
+        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
 
         if not adjoint:
             S_mDeriv = S_mDeriv(v)
@@ -309,6 +311,7 @@ class FieldsFDEM_h(FieldsFDEM):
         FieldsFDEM.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
+        self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
         self._MeMuI = self.survey.prob.MeMuI
         self._MfRho = self.survey.prob.MfRho
@@ -316,7 +319,7 @@ class FieldsFDEM_h(FieldsFDEM):
     def _hPrimary(self, hSolution, srcList):
         hPrimary = np.zeros_like(hSolution,dtype = complex)
         for i, src in enumerate(srcList):
-            hp = src.hPrimary(self.survey.prob)
+            hp = src.hPrimary(self.prob)
             if hp is not None:
                 hPrimary[:,i] += hp
             return hPrimary
@@ -337,7 +340,7 @@ class FieldsFDEM_h(FieldsFDEM):
     def _jPrimary(self, hSolution, srcList):
         jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]])
         for i, src in enumerate(srcList):
-            jp = src.jPrimary(self.survey.prob)
+            jp = src.jPrimary(self.prob)
             if jp is not None:
                 jPrimary[:,i] = jp 
         return jPrimary
@@ -345,7 +348,7 @@ class FieldsFDEM_h(FieldsFDEM):
     def _jSecondary(self, hSolution, srcList):
         j = self._edgeCurl*hSolution
         for i, src in enumerate(srcList):
-            _,S_e = src.eval(self.survey.prob)
+            _,S_e = src.eval(self.prob)
             if S_e is not None:
                 j[:,i] += -S_e
         return j
@@ -357,7 +360,7 @@ class FieldsFDEM_h(FieldsFDEM):
             return self._edgeCurl.T*v
 
     def _jSecondaryDeriv_m(self, src, v, adjoint=False):
-        _,S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
+        _,S_eDeriv = src.evalDeriv(self.prob, adjoint)
         S_eDeriv = S_eDeriv(v)
         if S_eDeriv is not None:
             return -S_eDeriv
diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py
index 44dbfa58..8a6937e1 100644
--- a/simpegEM/FDEM/SurveyFDEM.py
+++ b/simpegEM/FDEM/SurveyFDEM.py
@@ -1,7 +1,7 @@
 from SimPEG import Survey, Problem, Utils, np, sp
 from simpegEM.Utils import SrcUtils
 from simpegEM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
-
+from scipy.constants import mu_0
 
 ####################################################
 # Receivers 
@@ -189,11 +189,12 @@ class SrcFDEM_RawVec(SrcFDEM):
 class SrcFDEM_MagDipole(SrcFDEM):
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1.):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         self.moment = moment
+        self.mu = mu 
         SrcFDEM.__init__(self, rxList)
 
     def bPrimary(self,prob):
@@ -216,13 +217,13 @@ class SrcFDEM_MagDipole(SrcFDEM):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y')
+            a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
 
         else:
             srcfct = SrcUtils.MagneticDipoleVectorPotential
-            ax = srcfct(self.loc, gridX, 'x')
-            ay = srcfct(self.loc, gridY, 'y')
-            az = srcfct(self.loc, gridZ, 'z')
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
             a = np.concatenate((ax, ay, az))
 
         return C*a
@@ -235,17 +236,34 @@ class SrcFDEM_MagDipole(SrcFDEM):
         b_p = self.bPrimary(prob)
         return -1j*omega(self.freq)*b_p 
 
+    def S_e(self,prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
 
 
 class SrcFDEM_MagDipole_Bfield(SrcFDEM):
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
     #TODO: neither does moment
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1.):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         self.moment = moment
+        self.mu = mu
         SrcFDEM.__init__(self, rxList)
 
     def bPrimary(self,prob):
@@ -268,13 +286,13 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            bx = srcfct(self.loc, gridX, 'x')
-            bz = srcfct(self.loc, gridZ, 'z')
+            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
             b = np.concatenate((bx,bz))
         else:
-            bx = srcfct(self.loc, gridX, 'x')
-            by = srcfct(self.loc, gridY, 'y')
-            bz = srcfct(self.loc, gridZ, 'z')
+            bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
             b = np.concatenate((bx,by,bz))
 
         return b
@@ -287,14 +305,33 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
         b = self.bPrimary(prob)
         return -1j*omega(self.freq)*b
 
+    def S_e(self,prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
+
 
 class SrcFDEM_CircularLoop(SrcFDEM):
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
-    def __init__(self, rxList, freq, loc, orientation='Z', radius = 1.):
+    def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
         self.freq = float(freq)
         self.orientation = orientation
         self.radius = radius
+        self.mu = mu
+        self.loc = loc
         SrcFDEM.__init__(self, rxList)
 
     def bPrimary(self,prob):
@@ -316,25 +353,42 @@ class SrcFDEM_CircularLoop(SrcFDEM):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y', self.radius)
+            a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
             srcfct = SrcUtils.MagneticDipoleVectorPotential
-            ax = srcfct(self.loc, gridX, 'x', self.radius)
-            ay = srcfct(self.loc, gridY, 'y', self.radius)
-            az = srcfct(self.loc, gridZ, 'z', self.radius)
+            ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
+            ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
+            az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
             a = np.concatenate((ax, ay, az))
 
         return C*a
 
     def hPrimary(self,prob):
         b = self.bPrimary(prob)
-        return h_from_b
+        return 1./self.mu*b
 
     def S_m(self, prob):
         b = self.bPrimary(prob)
         return -1j*omega(self.freq)*b
 
+    def S_e(self,prob):
+        if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
+            return None
+        else:
+            eqLocs = prob._eqLocs
+
+            if eqLocs is 'FE':
+                mui_s = prob.curModel.mui - 1./self.mu
+                MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
+                C = prob.mesh.edgeCurl
+            elif eqLocs is 'EF':
+                mu_s = prob.curModel.mu - self.mu
+                MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
+                C = prob.mesh.edgeCurl.T
+
+            return -C.T * (MMui_s * self.bPrimary(prob))
+
 
 ####################################################
 # Survey
diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py
index e56a92ba..540ff24f 100644
--- a/simpegEM/Tests/test_FDEM.py
+++ b/simpegEM/Tests/test_FDEM.py
@@ -17,9 +17,11 @@ TOL = 1e-4
 FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
 CONDUCTIVITY = 1e1
 MU = mu_0
-freq = [1e-1, 2e-1]
+freq = 1e-1
 addrandoms = True
 
+SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+
 
 def getProblem(fdemType, comp):
     cs = 5.
@@ -35,23 +37,61 @@ def getProblem(fdemType, comp):
     x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
     XYZ = Utils.ndgrid(x,x,np.r_[0.])
     Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
-    Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.])
-    Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[1], loc=np.r_[0.,0.,0.])
-
-    survey = EM.FDEM.SurveyFDEM([Src0, Src1])
 
+    if SrcType is 'MagDipole':
+        Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
+    elif SrcType is 'MagDipole_Bfield':
+        Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
+    elif SrcType is 'CircularLoop':
+        Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
 
     if fdemType == 'e':
+        if SrcType is 'RawVec':
+            S_m = np.zeros(mesh.nF)
+            S_e = np.zeros(mesh.nE)
+            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
+
+        survey = EM.FDEM.SurveyFDEM([Src])
         prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
+
     elif fdemType == 'b':
+        if SrcType is 'RawVec':
+            S_m = np.zeros(mesh.nF)
+            S_e = np.zeros(mesh.nE)
+            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
+
+        survey = EM.FDEM.SurveyFDEM([Src])
         prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
+
     elif fdemType == 'j':
+        if SrcType is 'RawVec': 
+            S_m = np.zeros(mesh.nE)
+            S_e = np.zeros(mesh.nF)
+            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
+
+        survey = EM.FDEM.SurveyFDEM([Src])
         prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
+
     elif fdemType == 'h':
+        if SrcType is 'RawVec':
+            S_m = np.zeros(mesh.nE)
+            S_e = np.zeros(mesh.nF)
+            S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
+            S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
+            Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
+
+        survey = EM.FDEM.SurveyFDEM([Src])
         prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
+
     else:
         raise NotImplementedError()
     prb.pair(survey)
@@ -69,15 +109,15 @@ def adjointTest(fdemType, comp):
     print 'Adjoint %s formulation - %s' % (fdemType, comp)
 
     m  = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
-    mu = np.log(np.ones(prb.mesh.nC))*MU
+    mu = np.ones(prb.mesh.nC)*MU
 
     if addrandoms is True:
-        m  = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 
+        m  = m + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1 
         mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
 
-    prb.mu = mu 
     survey = prb.survey
-
+    prb.PropMap.PropModel.mu = mu
+    prb.PropMap.PropModel.mui = 1./mu
     u = prb.fields(m)
 
     v = np.random.rand(survey.nD)
@@ -99,10 +139,12 @@ def derivTest(fdemType, comp):
     mu = np.log(np.ones(prb.mesh.nC)*MU)
 
     if addrandoms is True:
-        x0  = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 
+        x0 = x0 + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1 
         mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
 
-    prb.mu = mu
+    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)
@@ -120,10 +162,11 @@ def crossCheckTest(fdemType, comp):
     mu = np.log(np.ones(mesh.nC)*MU)
 
     if addrandoms is True:
-        m  = m + np.random.randn(mesh.nC)*CONDUCTIVITY*1e-1 
+        m  = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 
         mu = mu + np.random.randn(mesh.nC)*MU*1e-1
 
-    prb1.mu = mu
+    prb1.PropMap.PropModel.mu = mu
+    prb1.PropMap.PropModel.mui = 1./mu
     survey1 = prb1.survey
     d1 = survey1.dpred(m)
 
diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py
index ec638a0f..1827f6b2 100644
--- a/simpegEM/Utils/SrcUtils.py
+++ b/simpegEM/Utils/SrcUtils.py
@@ -2,7 +2,7 @@ from SimPEG import *
 from scipy.special import ellipk, ellipe
 from scipy.constants import mu_0, pi
 
-def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
+def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMoment=(0., 0., 1.), mu = mu_0):
     """
         Calculate the vector potential of a set of magnetic dipoles
         at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
@@ -15,6 +15,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0
         :return: The vector potential each dipole at each observation location
     """
     #TODO: break this out!
+
     if type(component) in [list, tuple]:
         out = range(len(component))
         for i, comp in enumerate(component):
@@ -48,13 +49,13 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0
         dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0)
         mCr = np.cross(m, dR)
         r = np.sqrt((dR**2).sum(axis=1))
-        A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
+        A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3)
     if nSrc == 1:
         return A.flatten()
     return A
 
 
-def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
+def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0):
     """
         Calculate the vector potential of a set of magnetic dipoles
         at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
@@ -62,7 +63,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
         :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z)
         :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z)
         :param str component: The component to calculate - 'x', 'y', or 'z'
-        :param numpy.ndarray dipoleMoment: The vector dipole moment (vertical)
+        :param numpy.ndarray moment: The vector dipole moment (vertical)
         :rtype: numpy.ndarray
         :return: The vector potential each dipole at each observation location
     """
@@ -78,22 +79,22 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
 
     srcLoc = np.atleast_2d(srcLoc)
     obsLoc = np.atleast_2d(obsLoc)
-    dipoleMoment = np.atleast_2d(dipoleMoment)
+    moment = np.atleast_2d(moment)
 
     nFaces = obsLoc.shape[0]
     nSrc = srcLoc.shape[0]
 
-    m = np.array(dipoleMoment).repeat(nFaces, axis=0)
+    m = np.array(moment).repeat(nFaces, axis=0)
     B = np.empty((nFaces, nSrc))
     for i in range(nSrc):
         dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0)
         r = np.sqrt((dR**2).sum(axis=1))
         if dimInd == 0:
-            B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
+            B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
         elif dimInd == 1:
-            B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
+            B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
         elif dimInd == 2:
-            B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
+            B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
         else:
             raise Exception("Not Implemented")
     if nSrc == 1:
@@ -102,7 +103,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
 
 
 
-def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius):
+def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
     """
         Calculate the vector potential of horizontal circular loop
         at given locations
@@ -119,13 +120,13 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius):
     if type(component) in [list, tuple]:
         out = range(len(component))
         for i, comp in enumerate(component):
-            out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius)
+            out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
         return np.concatenate(out)
 
     if isinstance(obsLoc, Mesh.BaseMesh):
         mesh = obsLoc
         assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']"
-        return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius)
+        return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius, mu)
 
     srcLoc = np.atleast_2d(srcLoc)
     obsLoc = np.atleast_2d(obsLoc)