diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index 0eb5962f..c9020c80 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -29,7 +29,8 @@ class BaseFDEMProblem(BaseEMProblem):
             Ainv = self.Solver(A, **self.solverOpts)
             sol = Ainv * rhs
             Srcs = self.survey.getSrcByFreq(freq)
-            F[Srcs, self._fieldType] = sol
+            ftype = self._fieldType + '_sol'
+            F[Srcs, ftype] = sol
 
         return F
 
diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py
index 2e3e10d8..920cef67 100644
--- a/simpegEM/FDEM/FieldsFDEM.py
+++ b/simpegEM/FDEM/FieldsFDEM.py
@@ -9,10 +9,11 @@ class FieldsFDEM(Problem.Fields):
 
 
 class FieldsFDEM_e(FieldsFDEM):
-    knownFields = {'e':'E'}
+    knownFields = {'e_sol':'E'}
     aliasFields = {
-                    'b_sec' : ['e','F','_b_sec'],
-                    'b' : ['e','F','_b']
+                    'e' : ['e_sol','E','_e'],
+                    'b_sec' : ['e_sol','F','_b_sec'],
+                    'b' : ['e_sol','F','_b']
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -21,19 +22,31 @@ class FieldsFDEM_e(FieldsFDEM):
     def startup(self):
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
 
-    def _b_sec(self, e, src):
+    def _e(self, e_sol, src):
+        e = e_sol
+        e_p = src.e_p(self.survey.prob)
+        if e_p is not None:
+            e += e_p     
+        return e
+
+    def _b_sec(self, e_sol, src):
         C = self._edgeCurl
-        b_sec =  - 1./(1j*omega(src.freq))*(C * e)
+        b_sec =  - 1./(1j*omega(src.freq))*(C * e_sol)
         return b_sec
 
-    def _b_secDeriv(self, e, src, v, adjoint=False): 
+    def _b_secDeriv(self,e_sol, src, v, adjoint=False): 
         return None
 
-    def _b(self, e, src):
-        b = self._b_sec(e,src)
+    def _b(self, e_sol, src):
+        b = self._b_sec(e_sol, src)
         S_m, _ = src.eval(self.survey.prob)
         if S_m is not None:
             b += 1./(1j*omega(src.freq)) * S_m
+
+        b_p = src.b_p(self.survey.prob)
+        if b_p is not None:
+            b += b_p 
+
         return b
 
     def _bDeriv(self, e, src, v, adjoint=False):
@@ -50,10 +63,11 @@ class FieldsFDEM_e(FieldsFDEM):
 
 
 class FieldsFDEM_b(FieldsFDEM):
-    knownFields = {'b':'F'}
+    knownFields = {'b_sol':'F'}
     aliasFields = {
-                    'e_sec' : ['b','E','_e_sec'],
-                    'e' : ['b','E','_e']
+                    'b' : ['b_sol','F','_b'],
+                    'e_sec' : ['b_sol','E','_e_sec'],
+                    'e' : ['b_sol','E','_e']
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -66,22 +80,34 @@ class FieldsFDEM_b(FieldsFDEM):
         # self._getSource = self.survey.prob.getSource
         # self._getSourceDeriv = self.survey.prob.getSourceDeriv 
 
-    def _e_sec(self, b, src):
-        return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * b))
+    def _b(self, b_sol, src):
+        b = b_sol
+        b_p = src.b_p(self.survey.prob)
+        if b_p is not None:
+            b += b_p
+        return b  
 
-    def _e_secDeriv(self, b, src, v, adjoint=False):
+    def _e_sec(self, b_sol, src):
+        return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * b_sol))
+
+    def _e_secDeriv(self, b_sol, src, v, adjoint=False):
         return None
 
-    def _e(self, b, src):
-        e = self._e_sec(b,src)
+    def _e(self, b_sol, src):
+        e = self._e_sec(b_sol,src)
         _,S_e = src.eval(self.survey.prob)
         if S_e is not None:
-            e += S_e
+            e += -self._MeSigmaI*S_e
+
+        e_p = src.e_p(self.survey.prob)
+        if e_p is not None:
+            e += e_p
+
         return e
 
-    def _eDeriv(self, b, src, v, adjoint=False):
+    def _eDeriv(self, b_sol, src, v, adjoint=False):
         _,S_eDeriv = src.getSourceDeriv(self.survey.prob, v, adjoint)
-        e_secDeriv = self._e_secDeriv(b, src, v, adjoint)
+        e_secDeriv = self._e_secDeriv(b_sol, src, v, adjoint)
 
         if S_eDeriv is None & e_secDeriv is None:
             return None
@@ -94,10 +120,11 @@ class FieldsFDEM_b(FieldsFDEM):
 
 
 class FieldsFDEM_j(FieldsFDEM):
-    knownFields = {'j':'F'}
+    knownFields = {'j_sol':'F'}
     aliasFields = {
-                    'h_sec' : ['j','E','_h_sec'],
-                    'h' : ['j','E','_h']
+                    'j' : ['j_sol','F','_j'],
+                    'h_sec' : ['j_sol','E','_h_sec'],
+                    'h' : ['j_sol','E','_h']
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -111,10 +138,17 @@ class FieldsFDEM_j(FieldsFDEM):
         # self._getSourceDeriv = self.survey.prob.getSourceDeriv 
         self._curModel = self.survey.prob.curModel
 
-    def _h_sec(self, j, src): #v, adjoint=False
-        return - 1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfSigmai * j) ) 
+    def _j(self, j_sol, src):
+        j = j_sol
+        j_p = src.j_p(self.survey.prob) 
+        if j_p is not None:
+            j += j_p
+        return j
 
-    def _h_secDeriv(self, j, src, v, adjoint=False): 
+    def _h_sec(self, j_sol, src): #v, adjoint=False
+        return - 1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfSigmai * j_sol) ) 
+
+    def _h_secDeriv(self, j_sol, src, v, adjoint=False): 
         MeMuI = self._MeMuI
         C = self._edgeCurl
         sig = self._curModel.transform
@@ -128,16 +162,19 @@ class FieldsFDEM_j(FieldsFDEM):
         else:
             return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) )
 
-    def _h(self, j, src): #v, adjoint=False
-        h = self._h_sec(j,src)
+    def _h(self, j_sol, src): #v, adjoint=False
+        h = self._h_sec(j_sol,src)
         S_m,_ = src.eval(self.survey.prob)
         if S_m is not None:
             h += 1./(1j*omega(src.freq)) * self._MeMuI * S_m
+        h_p = src.h_p(self.survey.prob)
+        if h_p is not None:
+            h += h_p 
         return h
 
-    def _hDeriv(self, j, src, v, adjoint=False):
+    def _hDeriv(self, j_sol, src, v, adjoint=False):
         S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint)
-        h_secDeriv = self._h_secDeriv(j,src.freq, v, adjoint)
+        h_secDeriv = self._h_secDeriv(j_sol,src.freq, v, adjoint)
         if S_mDeriv is None & h_secDeriv is None:
             return None
         elif h_secDeriv is None:
@@ -147,11 +184,13 @@ class FieldsFDEM_j(FieldsFDEM):
         else:
             return 1./(1j*omega(src.freq)) * S_mDeriv + h_secDeriv
 
+
 class FieldsFDEM_h(FieldsFDEM):
-    knownFields = {'h':'E'}
+    knownFields = {'h_sol':'E'}
     aliasFields = {
-                    'j_sec' : ['h','F','_j_sec'],
-                    'j' : ['h','F','_j']
+                    'h' : ['h_sol','E','_h'],
+                    'j_sec' : ['h_sol','F','_j_sec'],
+                    'j' : ['h_sol','F','_j']
                   }
 
     def __init__(self,mesh,survey,**kwargs):
@@ -162,20 +201,30 @@ class FieldsFDEM_h(FieldsFDEM):
         self._MeMuI = self.survey.prob.MeMuI
         self._MfSigmai = self.survey.prob.MfSigmai
 
-    def _j_sec(self, h, src): # adjoint=False
-        return self._edgeCurl*h
+    def _h(self, h_sol, src):
+        h = h_sol
+        h_p = src.h_p(self.survey.prob)
+        if h_p is not None:
+            h += h_p
+        return h
 
-    def _j_secDeriv(self, h, src, v, adjoint=False): 
+    def _j_sec(self, h_sol, src): # adjoint=False
+        return self._edgeCurl*h_sol
+
+    def _j_secDeriv(self, h_sol, src, v, adjoint=False): 
         return None
 
-    def _j(self, h, src): # adjoint=False
-        j = self._j_sec(h,src)
+    def _j(self, h_sol, src): # adjoint=False
+        j = self._j_sec(h_sol,src)
         _,S_e = src.eval(self.survey.prob)
         if S_e is not None:
             j += -S_e
+        j_p = src.j_p(self.survey.prob)
+        if j_p is not None:
+            j += j_p 
         return j
 
-    def _jDeriv(self, h, src, v, adjoint=False):
+    def _jDeriv(self, h_sol, src, v, adjoint=False):
         _,S_eDeriv = src.getSourceDeriv(self.survey.prob, v, adjoint)
         j_secDeriv = self._j_secDeriv(j,src.freq, v, adjoint)
         if S_eDeriv is None & j_secDeriv is None:
diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py
index ddcf38d1..35f7ec97 100644
--- a/simpegEM/FDEM/SurveyFDEM.py
+++ b/simpegEM/FDEM/SurveyFDEM.py
@@ -1,6 +1,6 @@
 from SimPEG import Survey, Problem, Utils, np, sp
 from simpegEM.Utils import SrcUtils
-from simpegEM.Utils.EMUtils import omega
+from simpegEM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
 
 
 ####################################################
@@ -99,15 +99,58 @@ class SrcFDEM(Survey.BaseSrc):
         S_m = self._getS_m(prob)
         S_e = self._getS_e(prob)
 
-        if S_m is not None:
-            if len(S_m.shape) == 1: S_m = Utils.mkvc(S_m,2)
-        if S_e is not None:
-            if len(S_e.shape) == 1: S_e = Utils.mkvc(S_e,2)
+        if S_m is not None and S_m.ndim == 1: S_m = Utils.mkvc(S_m,2)
+        if S_e is not None and S_e.ndim == 1: S_e = Utils.mkvc(S_e,2)
+
         return S_m, S_e 
 
     def evalDeriv(self, prob, v, adjoint=None):
         return self._getS_mDeriv(prob,v,adjoint), self._getS_eDeriv(prob,v,adjoint)
 
+    def b_p(self,prob):
+        b_p = self._getb_p(prob)
+        if b_p is not None and b_p.ndim == 1: b_p = Utils.mkvc(b_p,2)
+        return b_p 
+
+    def h_p(self,prob):
+        h_p = self._geth_p(prob)
+        if h_p is not None and h_p.ndim == 1: h_p = Utils.mkvc(h_p,2)
+        return h_p
+
+    def e_p(self,prob):
+        e_p = self._gete_p(prob)
+        if e_p is not None and e_p.ndim == 1: e_p = Utils.mkvc(e_p,2)
+        return e_p
+
+    def j_p(self,prob):
+        j_p = self._getj_p(prob)
+        if j_p is not None and j_p.ndim == 1: j_p = Utils.mkvc(j_p,2)
+        return j_p
+
+    def _getb_p(self,prob):
+        return None
+
+    def _geth_p(self,prob):
+        return None
+
+    def _gete_p(self,prob):
+        return None
+
+    def _getj_p(self,prob):
+        return None
+
+    def _getS_m(self,prob):
+        return None
+
+    def _getS_e(self,prob):
+        return None
+
+    def _getS_mDeriv(self, prob, v, adjoint = False):
+        return None
+
+    def _getS_eDeriv(self, prob, v, adjoint = False):
+        return None
+
 
 class SrcFDEM_RawVec_e(SrcFDEM):
     """
@@ -123,18 +166,9 @@ class SrcFDEM_RawVec_e(SrcFDEM):
         self.freq = float(freq)
         SrcFDEM.__init__(self, rxList)
 
-    def _getS_m(self, prob):
-        return None 
-
     def _getS_e(self, prob):
         return self.S_e
 
-    def _getS_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def _getS_eDeriv(self, prob, v, adjoint = False):
-        return None
-
 
 class SrcFDEM_RawVec_m(SrcFDEM):
     """
@@ -153,15 +187,6 @@ class SrcFDEM_RawVec_m(SrcFDEM):
     def _getS_m(self, prob):
         return self.S_m
 
-    def _getS_e(self, prob):
-        return None
-
-    def _getS_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def _getS_eDeriv(self, prob, v, adjoint = False):
-        return None
-
 
 class SrcFDEM_RawVec(SrcFDEM):
     """
@@ -184,12 +209,7 @@ class SrcFDEM_RawVec(SrcFDEM):
     def _getS_e(self,prob):
         return self.S_e
 
-    def _getS_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def _getS_eDeriv(self, prob, v, adjoint = False):
-        return None
-
+ 
 class SrcFDEM_MagDipole(SrcFDEM):
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
@@ -200,7 +220,7 @@ class SrcFDEM_MagDipole(SrcFDEM):
         self.moment = moment
         SrcFDEM.__init__(self, rxList)
 
-    def _getS_m(self,prob):
+    def _getb_p(self,prob):
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -229,15 +249,16 @@ class SrcFDEM_MagDipole(SrcFDEM):
             az = srcfct(self.loc, gridZ, 'z')
             a = np.concatenate((ax, ay, az))
 
-        S_m = -1j*omega(self.freq)*C*a
+        return C*a
 
-        return S_m
+    def _geth_p(self,prob):
+        b = self._getb_p(prob)
+        return h_from_b(prob,b)
 
-    def _getS_e(self,prob):
-        return None
+    def _getS_m(self,prob):
+        b_p = self._getb_p(prob)
+        return -1j*omega(self.freq)*b_p 
 
-    def getSourceDeriv(self, prob, v, adjoint=None):
-        return None, None
 
 
 class SrcFDEM_MagDipole_Bfield(SrcFDEM):
@@ -251,7 +272,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
         self.moment = moment
         SrcFDEM.__init__(self, rxList)
 
-    def _getS_m(self,prob):
+    def _getb_p(self,prob):
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -280,17 +301,16 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM):
             bz = srcfct(self.loc, gridZ, 'z')
             b = np.concatenate((bx,by,bz))
 
+        return b
+
+    def _geth_p(self,prob):
+        b = self._getb_p(prob)
+        return h_from_b(prob, b)
+
+    def _getS_m(self,prob):
+        b = self._getb_p(prob)
         return -1j*omega(self.freq)*b
 
-    def _getS_e(self,prob):
-        return None
-
-    def _getS_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def _getS_eDeriv(self, prob, v, adjoint = False):
-        return None
-
 
 class SrcFDEM_CircularLoop(SrcFDEM):
 
@@ -301,14 +321,7 @@ class SrcFDEM_CircularLoop(SrcFDEM):
         self.radius = radius
         SrcFDEM.__init__(self, rxList)
 
-    def _getS_mDeriv(self, prob, v, adjoint = False):
-        return None
-
-    def _getS_eDeriv(self, prob, v, adjoint = False):
-        return None
-
-
-    def getSource(self, prob):
+    def _getb_p(self,prob):
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -336,7 +349,15 @@ class SrcFDEM_CircularLoop(SrcFDEM):
             az = srcfct(self.loc, gridZ, 'z', self.radius)
             a = np.concatenate((ax, ay, az))
 
-        return -1j*omega(self.freq)*C*a
+        return C*a
+
+    def _geth_p(self,prob):
+        b = self._getb_p(prob)
+        return h_from_b
+
+    def _getS_m(self, prob):
+        b = self._getb_p(prob)
+        return -1j*omega(self.freq)*b
 
 
 ####################################################
diff --git a/simpegEM/Utils/EMUtils.py b/simpegEM/Utils/EMUtils.py
index edd1d247..5cdaf150 100644
--- a/simpegEM/Utils/EMUtils.py
+++ b/simpegEM/Utils/EMUtils.py
@@ -1,10 +1,46 @@
 import numpy as np
 from scipy.constants import mu_0, epsilon_0
 
+# useful params
 def omega(freq):
     """Angular frequency, omega"""
     return 2.*np.pi*freq
 
 def k(freq, sigma, mu=mu_0, eps=epsilon_0):
-	w = omega(freq)
-	return np.sqrt(mu * eps * w**2 - 1j * w* mu * sigma)
\ No newline at end of file
+    w = omega(freq)
+    return np.sqrt(mu * eps * w**2 - 1j * w* mu * sigma)
+
+# Constitutive relations
+def e_from_j(prob,j):
+    eqLocs = prob._eqLocs
+    if eqLocs is 'FE':
+        MSigmaI = prob.MeSigmaI
+    elif eqLocs is 'EF':
+        MSigmaI = prob.MfSigmai
+    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.MfSigma
+    return MSigma*e
+
+def b_from_h(prob,h):
+    eqLocs = prob._eqLocs
+    if eqLocs is 'FE':
+        MMu = prob.MfMu
+    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 
+
+