diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py
index af2fb47b..1e53f12e 100644
--- a/simpegPF/BaseMag.py
+++ b/simpegPF/BaseMag.py
@@ -49,6 +49,9 @@ class BaseMagData(Data.BaseData):
         return bfz
 
         # return np.sqrt(bfx**2 + bfy**2 + bfz**2)
+    
+
+        # return np.sqrt(bfx**2 + bfy**2 + bfz**2)        
 
     @Utils.count
     def projectFieldsDeriv(self, B):
@@ -69,7 +72,7 @@ class BaseMagData(Data.BaseData):
         bfy = self.Qfy*B
         bfz = self.Qfz*B
 
-        return np.c_[bfx, bfy, bfz]
+        return np.r_[bfx, bfy, bfz]
 
 class MagDataBx(object):
     """docstring for MagDataBx"""
diff --git a/simpegPF/MagAnalytics.py b/simpegPF/MagAnalytics.py
index 2345d511..ec4d6485 100644
--- a/simpegPF/MagAnalytics.py
+++ b/simpegPF/MagAnalytics.py
@@ -120,7 +120,7 @@ def CongruousMagBC(mesh, Bo, chi):
 
     Bbcz = const/(rfun(mesh.gridFz[(indzd|indzu),:])**3)*(3*mdotrz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:])-mz)
 
-    return np.r_[Bbcx, Bbcy, Bbcz]
+    return np.r_[Bbcx, Bbcy, Bbcz], (1/gamma-1/(3+gamma))*1/V
 
 
 def MagSphereAnalFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
@@ -179,6 +179,15 @@ def MagSphereAnalFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
 
     return Bx, By, Bz
 
+def IDTtoxyz(Inc, Dec, Btot):
+    """
+        Convert from Inclination, Declination, Total  intensity of earth field to x, y, z
+    """
+    Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
+    By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
+    Bz = -Btot*np.sin(Inc/180.*np.pi)
+
+    return np.r_[Bx, By, Bz]
 
 if __name__ == '__main__':
 
@@ -195,7 +204,7 @@ if __name__ == '__main__':
     sph_ind = spheremodel(M3, 0, 0, 0, 100)
     chi[sph_ind] = chiblk
     mu = (1.+chi)*mu0
-    Bbc = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
+    Bbc, const = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
 
     flag = 'secondary'
     Box = 1.
diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py
index 9fa77b5d..eccf6786 100644
--- a/simpegPF/Magnetics.py
+++ b/simpegPF/Magnetics.py
@@ -6,7 +6,9 @@ from MagAnalytics import spheremodel, CongruousMagBC
 
 
 class MagneticsDiffSecondary(Problem.BaseProblem):
-    """Secondary field approach using differential equations!"""
+    """
+        Secondary field approach using differential equations!
+    """
 
     dataPair = BaseMag.BaseMagData
     modelPair = BaseMag.BaseMagModel
@@ -32,10 +34,12 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
 
     def makeMassMatrices(self, m):
         mu = self.model.transform(m)
-        self._MfMui = self.mesh.getFaceInnerProduct(1./mu)
+        self._MfMui = self.mesh.getFaceMass(1./mu)
+        # self._MfMui = self.mesh.getFaceInnerProduct(1./mu)
         #TODO: this will break if tensor mu
         self._MfMuI = Utils.sdiag(1./self._MfMui.diagonal())
-        self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0)
+        self._MfMu0 = self.mesh.getFaceMass(1/mu_0)
+        # self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0)
 
     def getB0(self):
         b0 = self.data.B0
@@ -51,10 +55,12 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
         Mc = Utils.sdiag(self.mesh.vol)
 
         chi = self.model.transform(m, asMu=False)
-        Bbc = CongruousMagBC(self.mesh, self.data.B0, chi)
-
+        Bbc, const = CongruousMagBC(self.mesh, self.data.B0, chi)
+        self.Bbc_const = const
+        self.Bbc = Bbc
         #TODO: put congrous BC back in
-        return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 #+ Mc*Dface*self._Pout.T*Bbc
+        # return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 #+ Mc*Dface*self._Pout.T*Bbc
+        return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc
 
     def getA(self, m):
         """
@@ -62,10 +68,13 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
 
         The A matrix has the form:
 
-        .. math::
+            .. math ::
 
-            \mathbf{A} = \mathbf{D}\mu\mathbf{G}u
+                \mathbf{A}\mathbf{u} = \mathbf{rhs}
 
+                \mathbf{A} = - \Div(\MfMui)^{-1}\Div^{T}
+                
+                \mathbf{rhs} = - \Div(\MfMui)^{-1}\mathbf{M}^f_{\\frac{1}{\mu_0}}\mathbf{B}_0 + \Div\mathbf{B}_0-\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}
 
 
         """
@@ -78,27 +87,70 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
         rhs = self.getRHS(m)
 
         m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/A.diagonal()))
-        phi, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1)
+        u, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1)
 
         B0 = self.getB0()
 
-        B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*phi
+        B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u
 
         #TODO: Create a mag fields object class.
         # F = self.getInitialFields()
-        # e.g. {'B': B, 'phi': phi}
-        return B
+        # e.g. {'B': B, 'u': u}
+
+        return {'B': B, 'u': u}
 
         # return self.forward(m, self.getRHS, self.calcFields, F=F)
 
     @Utils.timeIt
     def Jvec(self, m, v, u=None):
+        """
+            Computing Jacobian multiplied by vector
+            
+            By setting our problem as 
+
+            .. math ::
+
+                \mathbf{C}(\mathbf{m}, \mathbf{u}) = \mathbf{A}\mathbf{u} - \mathbf{rhs} = 0
+
+            And taking derivative w.r.t m
+
+            .. math ::
+
+                \\nabla \mathbf{C}(\mathbf{m}, \mathbf{u}) = \\nabla_m \mathbf{C}(\mathbf{m}) \delta \mathbf{m} +
+                                                             \\nabla_u \mathbf{C}(\mathbf{u}) \delta \mathbf{u} = 0
+
+                \\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = - [\\nabla_u \mathbf{C}(\mathbf{u})]^{-1}\\nabla_m \mathbf{C}(\mathbf{m})
+
+            With some linear algebra we can have 
+
+            .. math ::
+
+                \\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A}
+
+                \\nabla_m \mathbf{C}(\mathbf{m}) = 
+                \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}}                                
+
+            .. math :: 
+
+                \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} = 
+                \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[\Div \diag (\Div^T \mathbf{u}) \dMfMuI \\right]
+
+                \dMfMuI = \diag(\MfMui)^{-1}_{vec} \mathbf{Av}_{F2CC}^T\diag(\mathbf{v})\diag(\\frac{1}{\mu^2})
+
+                \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} =  \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[ 
+                \Div \diag(\M^f_{\mu_{0}^{-1} \mathbf{B}_0}) \dMfMuI \\right] - \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\\frac{\partial B_{sBC}}{\partial \mathbf{m}}
+
+
+        """
         if u is None:
             u = self.fields(m)
 
-        #TODO: B, phi = u['B'], u['phi']
+        #TODO: B, u = u['B'], u['u']
+
+        B, u = u['B'], u['u']
 
         mu = self.model.transform(m, asMu=True)
+        dmudm = self.model.transform(m, asMu=True)
 
         P = self.data.projectFieldsDeriv(u)
 
@@ -110,23 +162,30 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
 
         #TODO: only works for diagonal MfMui
         # Some chain rule!
-        harm_dm = Utils.sdiag(self.MfMui.diagonal()**(-2))
-        MfMu_dm = self.mesh.getFaceMassDeriv()
-        dmuI_dm = Utils.sdiag(mu**(-2))
-        mT_dm   = self.model.transformDeriv(m, asMu=True)
 
-        D = self._Div
+
+        # harm_dm = Utils.sdiag(self.MfMui.diagonal()**(-2))
+        # MfMu_dm = self.mesh.getFaceMassDeriv()
+        # dmuI_dm = Utils.sdiag(mu**(-2))
+        # mT_dm   = self.model.transformDeriv(m, asMu=True)
+
+        getFIPconst = 1./3
+        MfMuIvec = 1/self.MfMui.diagonal()*getFIPconst
+        dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(self.mesh.vol*1./mu**2)
+
+        Div = self._Div
         # lots-o-bracket for vector multiplication first!
-        MfMu_dmXv =  harm_dm * ( MfMu_dm * ( dmuI_dm * ( mT_dm * v ) ) )
-        dCdm_A = D * ( Utils.sdiag( D.T * phi ) * MfMu_dmXv )
+        # MfMu_dmXv =  harm_dm * ( MfMu_dm * ( dmuI_dm * ( mT_dm * v ) ) )
+        #dCdm_A = D * ( Utils.sdiag( D.T * u ) * MfMu_dmXv )
 
+        dCdm_A = dmudm*Div * ( Utils.sdiag( Div.T * u * dMfMuI ) )
 
         # rhs = D * MfMuI * MfMu0 * B0
 
         B0 = self.getB0()
 
         #TODO: add congrous stuff
-        dCdm_RHS = D * Utils.sdiag( self.MfMu0*B0  ) * MfMu_dmXv
+        dCdm_RHS = dmudm* Div * Utils.sdiag( self.MfMu0*B0  ) * dMfMuI - Utils.sdiag(self.mesh.vol)*self.mesh.faceDiv*self.Pout.T*self.Bbc*self.Bbc_const
 
 
         # c(m,u) = A(m)u - rhs(m)
@@ -134,9 +193,9 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
 
         solve = Solver(dCdu)
 
-        #TODO: Multiply by the dP(phi(m))/dphi
-        # We transformed phi in our fields object.
-        #         ( dBdphi *                   + dBdm(phi)  )
+        #TODO: Multiply by the dP(u(m))/du
+        # We transformed u in our fields object.
+        #         ( dBdu *                   + dBdm(u)  )
         Jv = - P *           solve.solve(dCdm)
         return Utils.mkvc(Jv)
 
@@ -177,24 +236,24 @@ if __name__ == '__main__':
     dpred = data.dpred(chi, u=B)
 
 
-    ##################
-    # Test J
-    ##################
+    # ##################
+    # # Test J
+    # ##################
 
-    d_chi = 0.8*chi #np.random.rand(mesh.nCz)
-    d_sph_ind = spheremodel(mesh, 0., 0., -100., 50)
-    d_chi[d_sph_ind] = 0.02
+    # d_chi = 0.8*chi #np.random.rand(mesh.nCz)
+    # d_sph_ind = spheremodel(mesh, 0., 0., -100., 50)
+    # d_chi[d_sph_ind] = 0.02
 
-    from SimPEG.Tests import checkDerivative
+    # from SimPEG.Tests import checkDerivative
 
-    derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)]
-    print '\n'
-    passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2)
+    # derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)]
+    # print '\n'
+    # passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2)
 
 
-    # plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))
+    # # plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))
 
-    # plt.show()
+    # # plt.show()