diff --git a/simpegEM/Base.py b/simpegEM/Base.py
index d1e6bb08..e88fdc5d 100644
--- a/simpegEM/Base.py
+++ b/simpegEM/Base.py
@@ -6,8 +6,6 @@ class BaseEMProblem(Problem.BaseProblem):
     def __init__(self, mesh, **kwargs):
         Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
-    solType = None
-    storeTheseFields = ['e', 'b']
 
     surveyPair = Survey.BaseSurvey
     dataPair = Survey.Data
@@ -105,5 +103,5 @@ class BaseEMProblem(Problem.BaseProblem):
 
     def fields(self, m):
         self.curModel = m
-        F = self.forward(m, self.getRHS, self.calcFields)
+        F = self.forward(m, self.getRHS)
         return F
diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index 95a66799..acd46d8e 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -1,56 +1,11 @@
 from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
-from SurveyFDEM import SurveyFDEM, FieldsFDEM
+from SurveyFDEM import SurveyFDEM
+from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b
 from simpegEM.Base import BaseEMProblem
 from simpegEM.Utils.EMUtils import omega
 
 
-
-# class FieldsTDEM_e_from_b(FieldsFDEM):
-#     """Fancy Field Storage for a TDEM survey."""
-#     knownFields = {'b_sec': 'F'}
-#     aliasFields = {
-#                     'b': ['b_sec','F','b_from_bsec'],
-#                     'e': ['b','E','e_from_b']
-#                   }
-
-#     def startup(self):
-#         self.MeSigmaI  = self.survey.prob.MeSigmaI
-#         self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
-#         self.MfMui     = self.survey.prob.MfMui
-
-#     def e_from_b(self, b, txInd, timeInd):
-#         # TODO: implement non-zero js
-#         return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
-
-#     def e_from_bDeriv(self, b, txInd, timeInd):
-#         # TODO: implement non-zero js
-#         return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
-
-
-#     def calcFields(self, sol, freq, fieldType, adjoint=False):
-#         e = sol
-#         if fieldType == 'e':
-#             return e
-#         elif fieldType == 'b':
-#             if not adjoint:
-#                 b = - self.mesh.edgeCurl * e 
-#                 b = 1./(1j*omega(freq)) * b
-#             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 BaseFDEMProblem(BaseEMProblem):
     """
         We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
@@ -62,26 +17,23 @@ class BaseFDEMProblem(BaseEMProblem):
 
     """
     surveyPair = SurveyFDEM
-    # fieldsPair = FieldsFDEM
+    fieldsPair = FieldsFDEM
 
-    def forward(self, m, RHS, CalcFields):
+    def forward(self, m, RHS):
 
-        # F = self.fieldsPair(self.mesh, self.survey)
-
-        F = FieldsFDEM(self.mesh, self.survey)
+        F = self.fieldsPair(self.mesh, self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
             rhs = RHS(freq)
             Ainv = self.Solver(A, **self.solverOpts)
             sol = Ainv * rhs
-            for fieldType in self.storeTheseFields:
-                Txs = self.survey.getTransmitters(freq)
-                F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
+            # for fieldType in self.storeTheseFields:
+            #     Txs = self.survey.getTransmitters(freq)
+            #     F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
 
-
-            # Txs = self.survey.getTransmitters(freq)
-            # F[Txs, 'e_sec'] = sol
+            Txs = self.survey.getTransmitters(freq)
+            F[Txs, self._fieldType] = sol
 
         return F
 
@@ -169,6 +121,10 @@ class BaseFDEMProblem(BaseEMProblem):
         return j_m, j_e
         # return np.concatenate(rhs).reshape((-1, len(Txs)), order='F') #, np.concatenate(j_e).reshape((-1, len(Txs)), order='F')
 
+    def getSourceDeriv(self,freq,adjoint=False):
+        return None, None
+
+
 ##########################################################################################
 ################################ E-B Formulation #########################################
 ##########################################################################################
@@ -191,9 +147,9 @@ class ProblemFDEM_e(BaseFDEMProblem):
 
 
     """
-    solType = 'e'
-    # _fieldType = 'e'
-    # fieldsPair = FieldsFDEM_e
+
+    _fieldType = 'e'
+    fieldsPair = FieldsFDEM_e
 
     def __init__(self, model, **kwargs):
         BaseFDEMProblem.__init__(self, model, **kwargs)
@@ -246,33 +202,34 @@ class ProblemFDEM_e(BaseFDEMProblem):
         return RHS
 
 
-    def calcFields(self, sol, freq, fieldType, adjoint=False):
-        e = sol
-        if fieldType == 'e':
-            return e
-        elif fieldType == 'b':
-            if not adjoint:
-                b = - self.mesh.edgeCurl * e 
-                b = 1./(1j*omega(freq)) * b
-            else:
-                b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
-            return b
-        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+    # def calcFields(self, sol, freq, fieldType, adjoint=False):
+    #     e = sol
+    #     if fieldType == 'e':
+    #         return e
+    #     elif fieldType == 'b':
+    #         if not adjoint:
+    #             b = - self.mesh.edgeCurl * e 
+    #             b = 1./(1j*omega(freq)) * b
+    #         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 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'
+    _fieldType = 'b'
+    fieldsPair = FieldsFDEM_b
 
     def __init__(self, model, **kwargs):
         BaseFDEMProblem.__init__(self, model, **kwargs)
@@ -345,40 +302,40 @@ class ProblemFDEM_b(BaseFDEMProblem):
         return RHS
 
 
-    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 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.curModel.transform
-            dsig_dm = self.curModel.transformDeriv
+    # def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
+    #     b = sol
+    #     if fieldType == 'e':
+    #         sig = self.curModel.transform
+    #         dsig_dm = self.curModel.transformDeriv
 
-            C = self.mesh.edgeCurl
-            mui = self.MfMui
+    #         C = self.mesh.edgeCurl
+    #         mui = self.MfMui
 
-            #TODO: This only works if diagonal (no tensors)...
-            dMeSigmaI_dI = - self.MeSigmaI**2
+    #         #TODO: This only works if diagonal (no tensors)...
+    #         dMeSigmaI_dI = - self.MeSigmaI**2
 
-            vec = C.T * ( mui * b )
-            dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(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)
+    #         vec = C.T * ( mui * b )
+    #         dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(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/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py
new file mode 100644
index 00000000..9c3c506f
--- /dev/null
+++ b/simpegEM/FDEM/FieldsFDEM.py
@@ -0,0 +1,120 @@
+from SimPEG import Survey, Problem, Utils, np, sp 
+from simpegEM.Utils.EMUtils import omega
+
+
+class FieldsFDEM(Problem.Fields):
+	"""Fancy Field Storage for a FDEM survey."""
+	knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'} # TODO: a, phi
+	dtype = complex
+
+	def calcFields(self,sol,txInd,freqInd,fieldType):
+		if fieldType == 'e':
+			return self._e(sol,txInd,freqInd)
+		elif fieldType == 'e_sec':
+			return self._e_sec(sol,txInd,freqInd)
+		elif fieldType == 'b':
+			return self._b(sol,txInd,freqInd)
+		elif fieldType == 'b_sec':
+			return self._b_sec(sol,txInd,freqInd)
+		else:
+			raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+	def calcFieldsDeriv(self,sol,txInd,freqInd,fieldType,adjoint=False):
+		if fieldType == 'e':
+			return self._eDeriv(sol,txInd,freqInd,adjoint)
+		elif fieldType == 'e_sec':
+			return self._e_secDeriv(sol,txInd,freqInd,adjoint)
+		elif fieldType == 'b':
+			return self._bDeriv(sol,txInd,freqInd,adjoint,adjoint)
+		elif fieldType == 'b_sec':
+			return self._b_secDeriv(sol,txInd,freqInd,adjoint,adjoint)
+		else:
+			raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+
+
+class FieldsFDEM_e(FieldsFDEM):
+	knownFields = {'e':'E'}
+	aliasFields = {
+					'b_sec' : ['e','F','_b_sec'],
+					'b' : ['e','F','_b']
+				  }
+
+	def __init__(self,mesh,survey,**kwargs):
+		FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+
+	def startup(self):
+		self.edgeCurl = self.survey.prob.mesh.edgeCurl
+		self.freqs = self.survey.freqs
+		self.getSource = self.survey.prob.getSource
+		self.getSourceDeriv = self.survey.prob.getSourceDeriv 
+
+	def _e(self, e, txInd, freqInd):
+		return e
+
+	def _eDeriv(self, e, txInd, freqInd, adjoint=False):
+		return None
+
+	def _b_sec(self, e, txInd, freqInd): #adjoint=False
+		iomegainv = 1./(1j*omega(self.freqs[freqInd]))
+		return -iomegainv * (self.edgeCurl * e)
+
+	def _b_secDeriv(self, e, txInd, freqInd, adjoint=False): 
+		return None
+
+	def _b(self, e, txInd, freqInd): #adjoint=False
+		freq = self.freqs[freqInd]
+		b_sec = self._bsec(e,txInd,freqInd)
+		j_m,_ = self.getSource(freq)
+		return 1./(1j*omega(freq)) + b_sec
+
+	def _bDreiv(self, e, txInd, freqInd, adjoint=False):
+		freq = self.freqs[freqInd]
+		j_mDeriv,_ = self.getSourceDeriv(freq, adjoint)
+		if j_mDeriv is None:
+			return None
+		else:
+			return 1./(1j*omega(freq)) * j_mDeriv
+
+
+class FieldsFDEM_b(FieldsFDEM):
+	knownFields = {'b':'F'}
+	aliasFields = {
+					'e_sec' : ['b','E','_e_sec'],
+					'e' : ['b','E','_e']
+				  }
+
+	def __init__(self,mesh,survey,**kwargs):
+		FieldsFDEM.__init__(self,mesh,survey,**kwargs)
+
+	def startup(self):
+		self.edgeCurl = self.survey.prob.mesh.edgeCurl
+		self.MeSigmaI = self.survey.prob.MeSigmaI
+		self.MfMui = self.survey.prob.MfMui
+		self.freqs = self.survey.freqs
+		self.getSource = self.survey.prob.getSource
+		self.getSourceDeriv = self.survey.prob.getSourceDeriv 
+
+	def _b(self, b, txInd, freqInd):
+		return b
+
+	def _bDeriv(self, b, txInd, freqInd, adjoint=False):
+		return None
+
+	def _e_sec(self, b, txInd, freqInd):
+		return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * b) )
+
+	def _e_secDeriv(self, b, txInd, freqInd, adjoint=False):
+		return None
+
+	def _e(self, b, txInd, freqInd):
+		e_sec = _e_sec(self, b, txInd, freqInd)
+		_, j_g = self.getSource(self.freqs[freqInd])
+		return e_s - j_g
+
+	def _eDeriv(self, b, txInd, freqInd, adjoint=False):
+		_,j_gDeriv = self.getSourceDeriv(self.freqs[freqInd], adjoint)
+		if j_gDeriv is None:
+			return None
+		else:
+			return -j_gDeriv	
\ No newline at end of file
diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py
index 67db17f0..d3fb9c90 100644
--- a/simpegEM/FDEM/SurveyFDEM.py
+++ b/simpegEM/FDEM/SurveyFDEM.py
@@ -2,9 +2,6 @@ from SimPEG import Survey, Problem, Utils, np, sp
 from simpegEM import Sources
 from simpegEM.Utils.EMUtils import omega
 
-def omega(freq):
-    """Change frequency to angular frequency, omega"""
-    return 2.*np.pi*freq
 
 class RxFDEM(Survey.BaseRx):
 
@@ -105,7 +102,7 @@ class TxFDEM(Survey.BaseTx):
 
         tx = self
         freq = tx.freq
-        solType = prob.solType
+        solType = prob._fieldType # Hack, should just ask whether j_m, j_g are defined on edges or faces
 
         if solType == 'e' or solType == 'b':
             gridEJx = prob.mesh.gridEx
@@ -212,12 +209,6 @@ class SimpleTxFDEM_m(TxFDEM):
         return self.vec, None
 
 
-class FieldsFDEM(Problem.Fields):
-    """Fancy Field Storage for a FDEM survey."""
-    knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'}
-    dtype = complex
-
-
 class SurveyFDEM(Survey.BaseSurvey):
     """
         docstring for SurveyFDEM
diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py
index 5892ce31..110b4d1e 100644
--- a/simpegEM/FDEM/__init__.py
+++ b/simpegEM/FDEM/__init__.py
@@ -1,2 +1,3 @@
 from SurveyFDEM import *
 from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
+from FieldsFDEM import *
\ No newline at end of file