diff --git a/SimPEG/EM/Analytics/FDEMcasing.py b/SimPEG/EM/Analytics/FDEMcasing.py
index b6c7acf1..d3db381d 100644
--- a/SimPEG/EM/Analytics/FDEMcasing.py
+++ b/SimPEG/EM/Analytics/FDEMcasing.py
@@ -1,6 +1,6 @@
 from SimPEG import Utils, np
 from scipy.constants import mu_0, epsilon_0
-from simpegEM.Utils.EMUtils import k
+from SimPEG.EM.Utils.EMUtils import k
 
 def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0):
     a = float(a)
@@ -16,43 +16,43 @@ def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps
 
     nobs = obsloc.shape[0]
     dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
-    
+
     r2 = _r2(dxyz[:,:2])
     sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
-    k2 = k(freq,sigma[2],mu[2],eps) 
-    
+    k2 = k(freq,sigma[2],mu[2],eps)
+
     return Kc1 * moment / (4.*np.pi) *np.exp(-1j*k2*sqrtr2z2) / sqrtr2z2
 
 
 def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
     HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
-    
+
     nobs = obsloc.shape[0]
     dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
-    
+
     r2 = _r2(dxyz[:,:2])
     sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
-    k2 = k(freq,sigma[2],mu[2],eps) 
-    
+    k2 = k(freq,sigma[2],mu[2],eps)
+
     return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
 
 
 def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
     HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
-    
+
     nobs = obsloc.shape[0]
     dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
-    
+
     r2z2 = _r2(dxyz)
     sqrtr2z2 = np.sqrt(r2z2)
-    k2 = k(freq,sigma[2],mu[2],eps) 
-    
+    k2 = k(freq,sigma[2],mu[2],eps)
+
     return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2)
 
 def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
     HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
     dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
-    
+
     nobs = obsloc.shape[0]
     dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
 
@@ -60,14 +60,14 @@ def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
     r = np.sqrt(r2)
     z = dxyz[:,2]
     sqrtr2z2 = np.sqrt(r2 + z**2)
-    k2 = k(freq,sigma[2],mu[2],eps) 
-    
+    k2 = k(freq,sigma[2],mu[2],eps)
+
     return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2)
 
 def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
     HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
     dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
-    
+
     nobs = obsloc.shape[0]
     dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
 
@@ -75,10 +75,10 @@ def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
     r = np.sqrt(r2)
     z = dxyz[:,2]
     sqrtr2z2 = np.sqrt(r2 + z**2)
-    k2 = k(freq,sigma[2],mu[2],eps) 
-    
+    k2 = k(freq,sigma[2],mu[2],eps)
+
     return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2)
-    
+
 def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
     return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
 
@@ -95,4 +95,4 @@ def getCasingBrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=eps
     return mu_0 * getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
 
 def getCasingBzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
-    return mu_0 * getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
\ No newline at end of file
+    return mu_0 * getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py
index c9296ec9..6d3df44f 100644
--- a/SimPEG/EM/FDEM/FDEM.py
+++ b/SimPEG/EM/FDEM/FDEM.py
@@ -2,33 +2,33 @@ from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyFDEM import SurveyFDEM
 from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j
-from simpegEM.Base import BaseEMProblem
-from simpegEM.Utils.EMUtils import omega
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.Utils.EMUtils import omega
 
 
 class BaseFDEMProblem(BaseEMProblem):
     """
-        We start by looking at Maxwell's equations in the electric 
-        field \\\(\\\mathbf{e}\\\) and the magnetic flux 
+        We start by looking at Maxwell's equations in the electric
+        field \\\(\\\mathbf{e}\\\) and the magnetic flux
         density \\\(\\\mathbf{b}\\\)
 
         .. math ::
 
             \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
             {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
-        
-        if using the E-B formulation (:code:`ProblemFDEM_e` 
-        or :code:`ProblemFDEM_b`) or the magnetic field 
+
+        if using the E-B formulation (:code:`ProblemFDEM_e`
+        or :code:`ProblemFDEM_b`) or the magnetic field
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
-        .. math :: 
+        .. math ::
 
             \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
         if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`).
 
-        The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) 
+        The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
 
     surveyPair = SurveyFDEM
@@ -36,7 +36,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
     def fields(self, m=None):
         """
-            Solve the forward problem for the fields. 
+            Solve the forward problem for the fields.
         """
 
         self.curModel = m
@@ -59,7 +59,7 @@ class BaseFDEMProblem(BaseEMProblem):
         """
 
         if f is None:
-           f = self.fields(m) 
+           f = self.fields(m)
 
         self.curModel = m
 
@@ -67,7 +67,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for freq in self.survey.freqs:
             dA_du = self.getA(freq) #
-            dA_duI = self.Solver(dA_du, **self.solverOpts) 
+            dA_duI = self.Solver(dA_du, **self.solverOpts)
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
@@ -97,7 +97,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         return Utils.mkvc(Jv)
 
-    def Jtvec(self, m, v, f=None): 
+    def Jtvec(self, m, v, f=None):
         """
             Sensitivity transpose times a vector
         """
@@ -157,7 +157,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
     def getSourceTerm(self, freq):
         """
-            Evaluates the sources for a given frequency and puts them in matrix form 
+            Evaluates the sources for a given frequency and puts them in matrix form
 
             :param float freq: Frequency
             :rtype: numpy.ndarray (nE or nF, nSrc)
@@ -165,11 +165,11 @@ class BaseFDEMProblem(BaseEMProblem):
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._eqLocs is 'FE':
-            S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) 
+            S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
             S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
         elif self._eqLocs is 'EF':
             S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
-            S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) 
+            S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
@@ -198,9 +198,9 @@ class ProblemFDEM_e(BaseFDEMProblem):
 
         .. math ::
 
-            \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}  
+            \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
 
-        which we solve for \\\(\\\mathbf{e}\\\). 
+        which we solve for \\\(\\\mathbf{e}\\\).
     """
 
     _fieldType = 'e'
@@ -271,10 +271,10 @@ class ProblemFDEM_e(BaseFDEMProblem):
                 return - 1j * omega(freq) * S_eDerivv
             else:
                 return None
-        else:   
+        else:
             S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
 
-            if S_mDerivv is not None and S_eDerivv is not None: 
+            if S_mDerivv is not None and S_eDerivv is not None:
                 return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv
             elif S_mDerivv is not None:
                 return C.T * (MfMui * S_mDerivv)
@@ -289,17 +289,17 @@ class ProblemFDEM_b(BaseFDEMProblem):
         We eliminate \\\(\\\mathbf{e}\\\) using
 
         .. math ::
-             
+
              \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right)
 
         and solve for \\\(\\\mathbf{b}\\\) using:
 
         .. math ::
-            
-            \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}  + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e} 
-    
+
+            \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}  + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
+
         .. note ::
-            The inverse problem will not work with full anisotropy 
+            The inverse problem will not work with full anisotropy
     """
 
     _fieldType = 'b'
@@ -345,14 +345,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
             return MeSigmaIDeriv.T * (C.T * v)
 
         if self._makeASymmetric is True:
-            return MfMui.T * ( C * ( MeSigmaIDeriv * v ) ) 
-        return C * ( MeSigmaIDeriv * v ) 
+            return MfMui.T * ( C * ( MeSigmaIDeriv * v ) )
+        return C * ( MeSigmaIDeriv * v )
 
 
     def getRHS(self, freq):
         """
             .. math ::
-                \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e} 
+                \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
 
             :param float freq: Frequency
             :rtype: numpy.ndarray (nE, nSrc)
@@ -404,7 +404,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
                 SrcDeriv = C * (self.MeSigmaI * S_eDeriv)
             elif adjoint:
                 SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv)
-        else: 
+        else:
             SrcDeriv = None
 
         if RHSderiv is not None and SrcDeriv is not None:
@@ -412,7 +412,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
         elif SrcDeriv is not None:
             RHSderiv = SrcDeriv
 
-        if RHSderiv is not None: 
+        if RHSderiv is not None:
             if self._makeASymmetric is True and not adjoint:
                 return MfMui.T * RHSderiv
 
@@ -427,13 +427,13 @@ class ProblemFDEM_b(BaseFDEMProblem):
 
 class ProblemFDEM_j(BaseFDEMProblem):
     """
-        We eliminate \\\(\\\mathbf{h}\\\) using 
+        We eliminate \\\(\\\mathbf{h}\\\) using
 
-        .. math :: 
+        .. math ::
 
             \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
 
-        and solve for \\\(\\\mathbf{j}\\\) using 
+        and solve for \\\(\\\mathbf{j}\\\) using
 
         .. math ::
 
@@ -510,7 +510,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
 
         S_m, S_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
-        MeMuI = self.MeMuI  
+        MeMuI = self.MeMuI
 
         RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
         if self._makeASymmetric is True:
@@ -521,7 +521,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
 
     def getRHSDeriv_m(self, src, v, adjoint=False):
         C = self.mesh.edgeCurl
-        MeMuI = self.MeMuI  
+        MeMuI = self.MeMuI
         S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
 
         if adjoint:
@@ -538,10 +538,10 @@ class ProblemFDEM_j(BaseFDEMProblem):
                 return - 1j * omega(freq) * S_eDerivv
             else:
                 return None
-        else:   
+        else:
             S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
 
-            if S_mDerivv is not None and S_eDerivv is not None: 
+            if S_mDerivv is not None and S_eDerivv is not None:
                 RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv
             elif S_mDerivv is not None:
                 RHSDeriv = C * (MeMuI * S_mDerivv)
@@ -560,17 +560,17 @@ class ProblemFDEM_j(BaseFDEMProblem):
 
 class ProblemFDEM_h(BaseFDEMProblem):
     """
-        We eliminate \\\(\\\mathbf{j}\\\) using 
+        We eliminate \\\(\\\mathbf{j}\\\) using
 
         .. math ::
 
             \mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e}
 
-        and solve for \\\(\\\mathbf{h}\\\) using 
+        and solve for \\\(\\\mathbf{h}\\\) using
 
         .. math ::
 
-            \\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} 
+            \\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
     """
 
@@ -612,7 +612,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
         """
             .. math ::
 
-                \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} 
+                \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
             :param float freq: Frequency
             :rtype: numpy.ndarray (nE, nSrc)
@@ -647,12 +647,12 @@ class ProblemFDEM_h(BaseFDEMProblem):
         S_eDeriv = S_eDeriv(v)
 
         if S_mDeriv is not None:
-            if RHSDeriv is not None: 
+            if RHSDeriv is not None:
                 RHSDeriv += S_mDeriv(v)
-            else: 
+            else:
                 RHSDeriv =  S_mDeriv(v)
         if S_eDeriv is not None:
-            if RHSDeriv is not None: 
+            if RHSDeriv is not None:
                 RHSDeriv += C.T * (MfRho * S_e)
             else:
                 RHSDeriv = C.T * (MfRho * S_e)
diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index ce655e2e..bb786bd1 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -1,5 +1,5 @@
-from SimPEG import Survey, Problem, Utils, np, sp 
-from simpegEM.Utils.EMUtils import omega
+from SimPEG import Survey, Problem, Utils, np, sp
+from SimPEG.EM.Utils.EMUtils import omega
 
 
 class FieldsFDEM(Problem.Fields):
@@ -30,11 +30,11 @@ class FieldsFDEM_e(FieldsFDEM):
         for i, src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
             if ep is not None:
-                ePrimary[:,i] = ep     
+                ePrimary[:,i] = ep
         return ePrimary
 
     def _eSecondary(self, eSolution, srcList):
-        return eSolution 
+        return eSolution
 
     def _e(self, eSolution, srcList):
         return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
@@ -54,7 +54,7 @@ class FieldsFDEM_e(FieldsFDEM):
                 bPrimary[:,i] += bp
         return bPrimary
 
-    def _bSecondary(self, eSolution, srcList): 
+    def _bSecondary(self, eSolution, srcList):
         C = self._edgeCurl
         b = (C * eSolution)
         for i, src in enumerate(srcList):
@@ -123,7 +123,7 @@ class FieldsFDEM_b(FieldsFDEM):
         return bSolution
 
     def _b(self, bSolution, srcList):
-        return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)  
+        return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
 
     def _bDeriv_u(self, src, v, adjoint=False):
         return None
@@ -142,7 +142,7 @@ class FieldsFDEM_b(FieldsFDEM):
 
     def _eSecondary(self, bSolution, srcList):
         e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
-        for i,src in enumerate(srcList): 
+        for i,src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
             if S_e is not None:
                 e[:,i] += -self._MeSigmaI * S_e
@@ -150,7 +150,7 @@ class FieldsFDEM_b(FieldsFDEM):
 
     def _eSecondaryDeriv_u(self, src, v, adjoint=False):
         if not adjoint:
-            return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) ) 
+            return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) )
         else:
             return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
 
@@ -168,7 +168,7 @@ class FieldsFDEM_b(FieldsFDEM):
 
         if not adjoint:
             de_dm = self._MeSigmaIDeriv(w) * v
-        elif adjoint: 
+        elif adjoint:
             de_dm = self._MeSigmaIDeriv(w).T * v
 
         _, S_eDeriv = src.evalDeriv(self.prob, adjoint)
@@ -205,7 +205,7 @@ class FieldsFDEM_j(FieldsFDEM):
         FieldsFDEM.__init__(self,mesh,survey,**kwargs)
 
     def startup(self):
-        self.prob = self.survey.prob 
+        self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
         self._MeMuI = self.survey.prob.MeMuI
         self._MfRho = self.survey.prob.MfRho
@@ -215,7 +215,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.prob) 
+            jp = src.jPrimary(self.prob)
             if jp is not None:
                 jPrimary[:,i] += jp
         return jPrimary
@@ -238,11 +238,11 @@ class FieldsFDEM_j(FieldsFDEM):
         for i, src in enumerate(srcList):
             hp = src.hPrimary(self.prob)
             if hp is not None:
-                hPrimary[:,i] = hp 
+                hPrimary[:,i] = hp
         return hPrimary
 
     def _hSecondary(self, jSolution, srcList):
-        h =  self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) ) 
+        h =  self._MeMuI * (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)
@@ -251,7 +251,7 @@ class FieldsFDEM_j(FieldsFDEM):
         return h
 
     def _hSecondaryDeriv_u(self, src, v, adjoint=False):
-        if not adjoint: 
+        if not adjoint:
             return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) )
         elif adjoint:
             return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v))
@@ -264,10 +264,10 @@ class FieldsFDEM_j(FieldsFDEM):
         MfRhoDeriv = self._MfRhoDeriv
         Me = self._Me
 
-        if not adjoint: 
-            hDeriv_m =  -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) 
+        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 ) ) 
+            hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
 
         S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
 
@@ -282,14 +282,14 @@ class FieldsFDEM_j(FieldsFDEM):
         return hDeriv_m
 
 
-    def _h(self, jSolution, srcList): 
+    def _h(self, jSolution, srcList):
         return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
 
     def _hDeriv_u(self, src, v, adjoint=False):
         return self._hSecondaryDeriv_u(src, v, adjoint)
 
     def _hDeriv_m(self, src, v, adjoint=False):
-        # assuming the primary doesn't depend on the model 
+        # assuming the primary doesn't depend on the model
         return self._hSecondaryDeriv_m(src, v, adjoint)
 
 
@@ -339,7 +339,7 @@ class FieldsFDEM_h(FieldsFDEM):
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
             if jp is not None:
-                jPrimary[:,i] = jp 
+                jPrimary[:,i] = jp
         return jPrimary
 
     def _jSecondary(self, hSolution, srcList):
@@ -370,5 +370,5 @@ class FieldsFDEM_h(FieldsFDEM):
         return self._jSecondaryDeriv_u(src,v,adjoint)
 
     def _jDeriv_m(self, src, v, adjoint=False):
-        # assuming the primary does not depend on the model 
-        return self._jSecondaryDeriv_m(src,v,adjoint)
\ No newline at end of file
+        # assuming the primary does not depend on the model
+        return self._jSecondaryDeriv_m(src,v,adjoint)
diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py
index 61711bdd..dbda9f80 100644
--- a/SimPEG/EM/FDEM/SurveyFDEM.py
+++ b/SimPEG/EM/FDEM/SurveyFDEM.py
@@ -1,10 +1,10 @@
 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 SimPEG.EM.Utils import SrcUtils
+from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
 from scipy.constants import mu_0
 
 ####################################################
-# Receivers 
+# Receivers
 ####################################################
 
 class RxFDEM(Survey.BaseRx):
@@ -99,13 +99,13 @@ class SrcFDEM(Survey.BaseSrc):
     def eval(self, prob):
         S_m = self.S_m(prob)
         S_e = self.S_e(prob)
-        return S_m, S_e 
+        return S_m, S_e
 
     def evalDeriv(self, prob, v, adjoint=False):
         return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
 
     def bPrimary(self, prob):
-        return None 
+        return None
 
     def hPrimary(self, prob):
         return None
@@ -225,7 +225,7 @@ class SrcFDEM_RawVec(SrcFDEM):
             return prob.Me * self._S_e
         return self._S_e
 
- 
+
 class SrcFDEM_MagDipole(SrcFDEM):
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
@@ -275,7 +275,7 @@ class SrcFDEM_MagDipole(SrcFDEM):
 
     def S_m(self, prob):
         b_p = self.bPrimary(prob)
-        return -1j*omega(self.freq)*b_p 
+        return -1j*omega(self.freq)*b_p
 
     def S_e(self, prob):
         if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py
index 8f4be4cd..5e82b7b2 100644
--- a/SimPEG/EM/TDEM/BaseTDEM.py
+++ b/SimPEG/EM/TDEM/BaseTDEM.py
@@ -1,10 +1,10 @@
 from SimPEG import Solver, Problem
 from SimPEG.Problem import BaseTimeProblem
-from simpegEM.Utils import SrcUtils
+from SimPEG.EM.Utils import SrcUtils
 from scipy.constants import mu_0
 from SimPEG.Utils import sdiag, mkvc
 from SimPEG import Utils, Mesh
-from simpegEM.Base import BaseEMProblem
+from SimPEG.EM.Base import BaseEMProblem
 import numpy as np
 
 
diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py
index bcd83962..9878438d 100644
--- a/SimPEG/EM/TDEM/SurveyTDEM.py
+++ b/SimPEG/EM/TDEM/SurveyTDEM.py
@@ -1,6 +1,6 @@
 from SimPEG import Utils, Survey, np
 from SimPEG.Survey import BaseSurvey
-from simpegEM.Utils import SrcUtils
+from SimPEG.EM.Utils import SrcUtils
 from BaseTDEM import FieldsTDEM