diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index 253b9984..8224f047 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -257,7 +257,7 @@ class Fields3D_e(Fields):
         """
 
         # assuming primary does not depend on the model
-        return Zero()
+        return src.ePrimaryDeriv(self.prob, v, adjoint) #Zero()
 
     def _bPrimary(self, eSolution, srcList):
         """
@@ -600,8 +600,8 @@ class Fields3D_b(Fields):
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
 
     def _j(self, bSolution, srcList):
         """
diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py
index 4590205d..b55ad1db 100644
--- a/SimPEG/EM/FDEM/SrcFDEM.py
+++ b/SimPEG/EM/FDEM/SrcFDEM.py
@@ -60,6 +60,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._bPrimary
 
+    def bPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary magnetic flux density
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def hPrimary(self, prob):
         """
         Primary magnetic field
@@ -72,6 +84,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._hPrimary
 
+    def hPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary magnetic field
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def ePrimary(self, prob):
         """
         Primary electric field
@@ -84,6 +108,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._ePrimary
 
+    def ePrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary electric field
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def jPrimary(self, prob):
         """
         Primary current density
@@ -96,6 +132,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._jPrimary
 
+    def jPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary current density
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def s_m(self, prob):
         """
         Magnetic source term
@@ -630,21 +678,118 @@ class PrimSecSigma(BaseSrc):
         return prob.MeSigmaDeriv(self.ePrimary(prob)) * v
 
 
-class PrimSecSigmaProblem(BaseSrc):
+class PrimSecMappedSigma(BaseSrc):
+
+    """
+    Primary-Secondary Source in which a mapping is provided to put the current model
+    onto the primary mesh. This is solved on every model update.
+
+    There are a lot of layers to the derivatives here!
+
+    **Required**
+    :param list rxList: Receiver List
+    :param float freq: frequency
+    :param ProblemFDEM primaryProblem: FDEM primary problem
+    :param SurveyFDEM primarySurvey: FDEM primary survey
+
+    **Optional**
+    :param Mapping map2meshs: mapping current model to act as primary model on the secondary mesh
+
+    """
 
     def __init__(self, rxList, freq, primaryProblem, primarySurvey, map2meshs = None ,**kwargs):
 
         self.primaryProblem = primaryProblem
         self.primarySurvey = primarySurvey
+
+        if self.primaryProblem.ispaired is False:
+            self.primaryProblem.pair(self.primarySurvey)
+
         self.map2meshs = map2meshs
 
         BaseSrc.__init__(self, rxList, freq=freq, **kwargs)
 
+    @property
+    def _ProjPrimary(self):
+        if getattr(self, '__ProjPrimary', None) is None:
+            self.__ProjPrimary = self.primaryProblem.mesh.getInterpolationMatCartMesh(meshs, locType='F', locTypeTo='E')
+        return self.__ProjPrimary
+
+
+    def _primaryFields(self, prob):
+
+        # TODO: cache and check if prob.curModel has changed
+        return self.primaryProblem.fields(prob.curModel)
+
+    def _primaryFieldsDeriv(self, prob, v, adjoint=False):
+        if adjoint:
+            raise NotImplementedError
+
+        # TODO: this should not be hard-coded for j
+        jp = self._primaryFields(prob)[:,'j']
+
+        # TODO: pull apart Jvec so that don't have to copy paste this code in
+        A = self.primaryProblem.getA(self.freq)
+        Ainv = self.primaryProblem.Solver(A, **self.primaryProblem.solverOpts) # create the concept of Ainv (actually a solve)
+
+        df_dm_v = np.zeros(self.primaryProblem.survey.nSrc,len(jp))
+
+        # TODO: this will probably break if we have more than one source
+        for i, src in enumerate(self.primaryProblem.survey.getSrcByFreq(freq)):
+            u_src = f[src, self.primaryProblem._solutionType]
+            dA_dm_v = self.primaryProblem.getADeriv(freq, u_src, v)
+            dRHS_dm_v = self.primaryProblem.getRHSDeriv(freq, src, v)
+            du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+            df_dmFun = getattr(f, '_{0}Deriv'.format('j'), None)
+            df_dm_v[i,:] += df_dmFun(src, du_dm_v, v, adjoint=False)
+            # Jv[src, :] = rx.evalDeriv(src, self.primaryProblem.mesh, f, df_dm_v)
+
+        Ainv.clean()
+
+        return Utils.mkvc(df_dm_v)
 
     def ePrimary(self, prob):
 
+        jp = self._primaryFields(prob)[:,'j']
+        ep = self.primaryProblem.MfI * (self.primaryProblem.MfRho * jp)
+        ep = self._ProjPrimary * ep
+
+        return ep
+
     def ePrimaryDeriv(self, prob, v, adjoint=False):
 
+        if adjoint:
+            raise NotImplementedError
+
+        jp = self._primaryFields(prob)[:,'j']
+
+        epDeriv = self._ProjPrimary * (
+                    self.primaryProblem.MfI * ( self.primaryProblem.MfRhoDeriv(jp) * v )
+                    +
+                    self._primaryFieldsDeriv(prob, v, adjoint)
+                    )
+
+        return epDeriv
+
+
+    def s_e(self, prob):
+        sigmaPrimary = self.map2meshs * prob.curModel
+        return (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimary(prob)
+
+
+    def s_eDeriv(self, prob, v, adjoint=False):
+        if adjoint:
+            raise NotImplementedError
+            return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
+
+        sigmaPrimary = self.map2meshs * prob.curModel
+        sigmaPrimaryDeriv = self.map2meshs.deriv(prob.curModel)
+        return (prob.MeSigmaDeriv(self.ePrimary(prob)) * v
+                - prob.mesh.getEdgeInnerProductDeriv(sigmaPrimary)(self.ePrimary(prob)) * sigmaPrimaryDeriv * v
+                + (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimaryDeriv(prob, v)
+                )
+