diff --git a/simpegMT/Base.py b/simpegMT/Base.py
deleted file mode 100644
index 6215f5f5..00000000
--- a/simpegMT/Base.py
+++ /dev/null
@@ -1,68 +0,0 @@
-from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
-from scipy.constants import mu_0
-
-class BaseMTProblem(Problem.BaseProblem):
-
-    def __init__(self, mesh, **kwargs):
-        Problem.BaseProblem.__init__(self, mesh, **kwargs)
-
-    solType = None
-    storeTheseFields = ['e', 'b']
-
-    surveyPair = Survey.BaseSurvey
-    dataPair = Survey.Data
-
-    Solver = SimpegSolver
-    solverOpts = {}
-
-    ####################################################
-    # Mass Matrices
-    ####################################################
-
-    @property
-    def MfMui(self):
-        #TODO: assuming constant mu
-        if getattr(self, '_MfMui', None) is None:
-            self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
-        return self._MfMui
-
-    @property
-    def Me(self):
-        if getattr(self, '_Me', None) is None:
-            self._Me = self.mesh.getEdgeInnerProduct()
-        return self._Me
-
-    @property
-    def MeSigma(self):
-        #TODO: hardcoded to sigma as the model
-        if getattr(self, '_MeSigma', None) is None:
-            sigma = self.curTModel
-            self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
-        return self._MeSigma
-
-    @property
-    def MeSigmaI(self):
-        #TODO: hardcoded to sigma as the model
-        if getattr(self, '_MeSigmaI', None) is None:
-            sigma = self.curTModel
-            self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
-        return self._MeSigmaI
-
-    curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.')
-
-    @property
-    def curTModel(self):
-        if getattr(self, '_curTModel', None) is None:
-            self._curTModel = self.mapping.transform(self.curModel)
-        return self._curTModel
-
-    @property
-    def curTModelDeriv(self):
-        if getattr(self, '_curTModelDeriv', None) is None:
-            self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
-        return self._curTModelDeriv
-
-    def fields(self, m):
-        self.curModel = m
-        F = self.forward(m, self.getRHS, self.calcFields)
-        return F
diff --git a/simpegMT/FDEM/__init__.py b/simpegMT/FDEM/__init__.py
deleted file mode 100644
index 70124686..00000000
--- a/simpegMT/FDEM/__init__.py
+++ /dev/null
@@ -1,2 +0,0 @@
-from SurveyFDEM import *
-from FDEM import ProblemFDEM_e, ProblemFDEM_b
diff --git a/simpegMT/FDEM/FDEM.py b/simpegMT/ProblemMT.py
similarity index 51%
rename from simpegMT/FDEM/FDEM.py
rename to simpegMT/ProblemMT.py
index e3886280..a6669fa7 100644
--- a/simpegMT/FDEM/FDEM.py
+++ b/simpegMT/ProblemMT.py
@@ -1,41 +1,213 @@
 from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
-from SurveyFDEM import SurveyFDEM, FieldsFDEM
-# from simpegMT.Utils import Sources
-from simpegMT.Base import BaseMTProblem
+from SurveyMT import SurveyMT, FieldsMT
+
 
 def omega(freq):
     """Change frequency to angular frequency, omega"""
     return 2.*np.pi*freq
 
-class BaseFDEMProblem(BaseMTProblem):
-    """
-        We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
 
-        .. math::
+class MTProblem(Problem.BaseProblem):
 
-            \\nabla \\times \\vec{E} + i \\omega \\vec{B} = 0 \\\\
-            \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s}
+    def __init__(self, mesh, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
-    """
+    solType = 'e'
+    storeTheseFields = ['e', 'b']
 
-    surveyPair = SurveyFDEM
+    surveyPair = SurveyMT
+    dataPair = Survey.Data
 
-    def forward(self, m, RHS, CalcFields):
+    Solver = SimpegSolver
+    solverOpts = {}
 
-        F = FieldsFDEM(self.mesh, self.survey)
+    ####################################################
+    # Mass Matrices
+    ####################################################
+
+    @property
+    def MfMui(self):
+        #TODO: assuming constant mu
+        if getattr(self, '_MfMui', None) is None:
+            self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
+        return self._MfMui
+
+    @property
+    def Me(self):
+        if getattr(self, '_Me', None) is None:
+            self._Me = self.mesh.getEdgeInnerProduct()
+        return self._Me
+
+    @property
+    def MeSigma(self):
+        #TODO: hardcoded to sigma as the model
+        if getattr(self, '_MeSigma', None) is None:
+            sigma = self.curTModel
+            self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
+        return self._MeSigma
+
+    # TODO:
+    # MeSigmaBG
+
+    @property
+    def MeSigmaI(self):
+        #TODO: hardcoded to sigma as the model
+        if getattr(self, '_MeSigmaI', None) is None:
+            sigma = self.curTModel
+            self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
+        return self._MeSigmaI
+
+    curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.')
+
+    @property
+    def curTModel(self):
+        if getattr(self, '_curTModel', None) is None:
+            self._curTModel = self.mapping.transform(self.curModel)
+        return self._curTModel
+
+    @property
+    def curTModelDeriv(self):
+        if getattr(self, '_curTModelDeriv', None) is None:
+            self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
+        return self._curTModelDeriv
+
+    def fields(self, m):
+        self.curModel = m
+        RHS, CalcFields = self.getRHS, self.calcFields
+
+        F = FieldsMT(self.mesh, self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
             rhs = RHS(freq)
-            solver = self.Solver(A, **self.solverOpts)
-            sol = solver.solve(rhs)
-            for fieldType in self.storeTheseFields:
-                Txs = self.survey.getTransmitters(freq)
-                F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e = Ainv * rhs # is this e?
+
+            Src = self.survey.getSources(freq)
+            F[Src, 'e'] = e
+            F[Src, 'b'] = self.mesh.edgeCurl * e # ???
 
         return F
 
+
+    def getA(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+        mui = self.MfMui
+        sig = self.MeSigma
+        C = self.mesh.edgeCurl
+
+        return C.T*mui*C + 1j*omega(freq)*sig
+    def getAbg(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+        mui = self.MfMui
+        sigBG = self.MeSigmaBG
+        C = self.mesh.edgeCurl
+
+        return C.T*mui*C + 1j*omega(freq)*sigBG
+
+    def getADeriv(self, freq, u, v, adjoint=False):
+        sig = self.curTModel
+        dsig_dm = self.curTModelDeriv
+        dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
+
+        if adjoint:
+            return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
+
+        return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
+
+    def getRHS(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2)
+            :return: one RHS for both polarizations
+        """
+        raise NotImplementedError()
+
+        getAbg(freq)
+
+
+        """
+        Put this in MT.Sources.EldadsSource
+
+        from simpegMT.Utils import get1DEfields
+        # Get a 1d solution for a halfspace background
+        mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
+        e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
+        # Setup x (east) polarization (_x)
+        ex_x = np.zeros(M.vnEx,dtype=complex)
+        ey_x = np.zeros((M.nEy,1),dtype=complex)
+        ez_x = np.zeros((M.nEz,1),dtype=complex)
+        # Assign the source to ex_x
+        for i in arange(M.vnEx[0]):
+            for j in arange(M.vnEx[2]):
+                ex_x[i,j,:] = e0_1d
+        eBG_x = np.vstack((simpeg.Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))
+        rhs_x = ABG.dot(eBG_x)
+        """
+
+        Txs = self.survey.getTransmitters(freq)
+
+        # assert that only one Tx/src?
+        # Create the two polarizations at this freq and return np array (nE,2).
+
+        # solve analytic.... get p1 p2
+
+        # Abg * [p1,p2] = rhs
+
+        rhs = range(len(Txs))
+        for i, tx in enumerate(Txs):
+            if tx.txType == 'VMD': # EH source.
+                src = Sources.MagneticDipoleVectorPotential # this is where you would put multiple types of boundary conditions.
+            else:
+                raise NotImplemented('%s txType is not implemented' % tx.txType)
+            SRCx = src(tx.loc, self.mesh.gridEx, 'x')
+            SRCy = src(tx.loc, self.mesh.gridEy, 'y')
+            SRCz = src(tx.loc, self.mesh.gridEz, 'z')
+            rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
+
+        a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
+        mui = self.MfMui
+        C = self.mesh.edgeCurl
+
+        j_s = C.T*mui*C*a
+        return -1j*omega(freq)*j_s
+
+    ##################################################################
+    # Inversion stuff
+    ##################################################################
+    # Not really used now....
+
+
+    def calcFields(self, sol, freq, fieldType, adjoint=False):
+        e = sol
+        if fieldType == 'e':
+            return e
+        elif fieldType == 'b':
+            if not adjoint:
+                b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl * e )
+            else:
+                b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
+            return b
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+    def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
+        e = sol
+        if fieldType == 'e':
+            return None
+        elif fieldType == 'b':
+            return None
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+
     def Jvec(self, m, v, u=None):
         if u is None:
            u = self.fields(m)
@@ -104,193 +276,3 @@ class BaseFDEMProblem(BaseMTProblem):
                         raise Exception('Must be real or imag')
 
         return Jtv
-
-
-class ProblemFDEM_e(BaseFDEMProblem):
-    """
-        By eliminating the magnetic flux density using
-
-        .. math::
-
-            \\vec{B} = \\frac{-1}{i\\omega}\\nabla\\times\\vec{E},
-
-        we can write Maxwell's equations as a second order system in \\ \\vec{E} \\ only:
-
-        .. math::
-
-            \\nabla \\times \\mu^{-1} \\nabla \\times \\vec{E} + i \\omega \\sigma \\vec{E} = \\vec{J_s}
-
-        This is the definition of the Forward Problem using the E-formulation of Maxwell's equations.
-
-
-    """
-    solType = 'e'
-
-    def __init__(self, model, **kwargs):
-        BaseFDEMProblem.__init__(self, model, **kwargs)
-
-    def getA(self, freq):
-        """
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        mui = self.MfMui
-        sig = self.MeSigma
-        C = self.mesh.edgeCurl
-
-        return C.T*mui*C + 1j*omega(freq)*sig
-
-    def getADeriv(self, freq, u, v, adjoint=False):
-        sig = self.curTModel
-        dsig_dm = self.curTModelDeriv
-        dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
-
-        if adjoint:
-            return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
-
-        return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
-
-    def getRHS(self, freq):
-        """
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nTx)
-            :return: RHS
-        """
-        Txs = self.survey.getTransmitters(freq)
-        rhs = range(len(Txs))
-        for i, tx in enumerate(Txs):
-            if tx.txType == 'VMD':
-                src = Sources.MagneticDipoleVectorPotential
-            else:
-                raise NotImplemented('%s txType is not implemented' % tx.txType)
-            SRCx = src(tx.loc, self.mesh.gridEx, 'x')
-            SRCy = src(tx.loc, self.mesh.gridEy, 'y')
-            SRCz = src(tx.loc, self.mesh.gridEz, 'z')
-            rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
-
-        a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
-        mui = self.MfMui
-        C = self.mesh.edgeCurl
-
-        j_s = C.T*mui*C*a
-        return -1j*omega(freq)*j_s
-
-    def calcFields(self, sol, freq, fieldType, adjoint=False):
-        e = sol
-        if fieldType == 'e':
-            return e
-        elif fieldType == 'b':
-            if not adjoint:
-                b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl * e )
-            else:
-                b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
-            return b
-        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
-
-    def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
-        e = sol
-        if fieldType == 'e':
-            return None
-        elif fieldType == 'b':
-            return None
-        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
-
-
-class ProblemFDEM_b(BaseFDEMProblem):
-    """
-        Solving for b!
-    """
-    solType = 'b'
-
-    def __init__(self, model, **kwargs):
-        BaseFDEMProblem.__init__(self, model, **kwargs)
-
-    def getA(self, freq):
-        """
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        mui = self.MfMui
-        sigI = self.MeSigmaI
-        C = self.mesh.edgeCurl
-
-        return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui
-
-    def getADeriv(self, freq, u, v, adjoint=False):
-
-        mui = self.MfMui
-        C = self.mesh.edgeCurl
-        sig = self.curTModel
-        dsig_dm = self.curTModelDeriv
-        #TODO: This only works if diagonal (no tensors)...
-        dMeSigmaI_dI = - self.MeSigmaI**2
-
-        vec = (C.T*(mui*u))
-        dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=vec)
-
-        if adjoint:
-            return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) )
-
-        return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
-
-    def getRHS(self, freq):
-        """
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nTx)
-            :return: RHS
-        """
-        Txs = self.survey.getTransmitters(freq)
-        rhs = range(len(Txs))
-        for i, tx in enumerate(Txs):
-            if tx.txType == 'VMD':
-                src = Sources.MagneticDipoleVectorPotential
-            else:
-                raise NotImplemented('%s txType is not implemented' % tx.txType)
-            SRCx = src(tx.loc, self.mesh.gridEx, 'x')
-            SRCy = src(tx.loc, self.mesh.gridEy, 'y')
-            SRCz = src(tx.loc, self.mesh.gridEz, 'z')
-            rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
-
-        a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
-        mui = self.MfMui
-        C = self.mesh.edgeCurl
-
-        b_0 = C*a
-        return -1j*omega(freq)*mui*b_0
-
-    def calcFields(self, sol, freq, fieldType, adjoint=False):
-        b = sol
-        if fieldType == 'e':
-            if not adjoint:
-                e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) )
-            else:
-                e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) )
-            return e
-        elif fieldType == 'b':
-            return b
-        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
-
-    def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
-        b = sol
-        if fieldType == 'e':
-            sig = self.curTModel
-            dsig_dm = self.curTModelDeriv
-
-            C = self.mesh.edgeCurl
-            mui = self.MfMui
-
-            #TODO: This only works if diagonal (no tensors)...
-            dMeSigmaI_dI = - self.MeSigmaI**2
-
-            vec = C.T * ( mui * b )
-            dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=vec)
-            if not adjoint:
-                return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) )
-            else:
-                return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) )
-        elif fieldType == 'b':
-            return None
-        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
-
diff --git a/simpegMT/FDEM/SurveyFDEM.py b/simpegMT/SurveyMT.py
similarity index 92%
rename from simpegMT/FDEM/SurveyFDEM.py
rename to simpegMT/SurveyMT.py
index 67db6674..4072a4ad 100644
--- a/simpegMT/FDEM/SurveyFDEM.py
+++ b/simpegMT/SurveyMT.py
@@ -65,17 +65,19 @@ class RxFDEM(Survey.BaseRx):
         return Pv
 
 
+# Call this Source or polarization or something...?
 class TxFDEM(Survey.BaseTx):
 
     freq = None #: Frequency (float)
 
     rxPair = RxFDEM
 
-    knownTxTypes = ['VMD']
+    knownTxTypes = ['VMD'] # Polarization...
 
-    def __init__(self, loc, txType, freq, rxList):
+    def __init__(self, loc, txType, freq, rxList): # remove txType? hardcode to one thing. always polarizations
         self.freq = float(freq)
         Survey.BaseTx.__init__(self, loc, txType, rxList)
+        # Survey.BaseTx.__init__(self, loc, 'polarization', rxList)
 
 
 
@@ -124,6 +126,7 @@ class SurveyFDEM(Survey.BaseSurvey):
                 self._nTxByFreq[freq] = len(self.getTransmitters(freq))
         return self._nTxByFreq
 
+    # TODO: Rename to getSources
     def getTransmitters(self, freq):
         """Returns the transmitters associated with a specific frequency."""
         assert freq in self._freqDict, "The requested frequency is not in this survey."