removed MT1d example from this pr

if mumps import fails use LU

SimPEG/DataMisfit.py

@@ -59,20 +59,6 @@ class BaseDataMisfit(object):
         """
         raise NotImplementedError('This method should be overwritten.')
 
-    # TODO: implement target misfit as a property, or possibly as an inversion directive.
-
-    # def target(self, forward):
-    #     """target(forward)
-
-    #         Target for data misfit. By default this is the number of data,
-    #         which satisfies the Discrepancy Principle.
-
-    #         :rtype: float
-    #         :return: data misfit target
-
-    #     """
-    #     prob, survey = self.splitForward(forward)
-    #     return survey.nD
 
 
 class l2_DataMisfit(BaseDataMisfit):
@@ -103,10 +89,18 @@ class l2_DataMisfit(BaseDataMisfit):
         """
 
         if getattr(self, '_Wd', None) is None:
-            print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.'
+            
             survey = self.survey
-            eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
-            self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
+
+            if getattr(survey,'std', None) is None:
+                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
+                survey.std = 0.05
+
+            if getattr(survey, 'eps', None) is None:
+                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
+                survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
+
+            self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
         return self._Wd
 
     @Wd.setter

SimPEG/Directives.py

@@ -239,51 +239,14 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 
 
 
-class update_IRLS(InversionDirective):
+# class UpdateReferenceModel(Parameter):
 
-     m = None
-     eps_min = None
-     factor = None
-     gamma = None
-     phi_m_last = None
-     
-     def initialize(self):
-        
-         # Scale the regularization for changes in norm
-         if getattr(self, 'phi_m_last', None) is not None:
-             self.reg.gamma = 1.
-             phim_new = self.reg.eval(self.invProb.curModel)
-             self.gamma = self.phi_m_last / phim_new
-         
-         self.reg.gamma = self.gamma
-         
-     def endIter(self):
-         # Cool the threshold parameter 
-         if getattr(self, 'factor', None) is not None:
-             eps = self.reg.eps / self.factor
-             
-             if getattr(self, 'eps_min', None) is not None:
-                 self.reg.eps = np.max([self.eps_min,eps])
-             else:
-                 self.reg.eps = eps
-                 
-         
-         # Update the model used for the IRLS weights
-         if getattr(self, 'm', None) is None:
-             self.reg.m = self.invProb.curModel
-         
-         # Update the pre-conditioner
-         diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.prob.mesh.nC))**2.
-         PC     = Utils.sdiag(diagA**-1.)
+#     mref0 = None
 
-         self.opt.approxHinv = PC
-         
-         phim_new = self.reg.eval(self.invProb.curModel)
-         self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new
-
-#==============================================================================
-#          import pylab as plt
-#          plt.figure()
-#          ax = plt.subplot(221)
-#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
-#==============================================================================
+#     def nextIter(self):
+#         mref = getattr(self, 'm_prev', None)
+#         if mref is None:
+#             if self.debug: print 'UpdateReferenceModel is using mref0'
+#             mref = self.mref0
+#         self.m_prev = self.invProb.m_current
+#         return mref

SimPEG/EM/FDEM/FDEM.py

@@ -15,18 +15,20 @@ class BaseFDEMProblem(BaseEMProblem):
         .. math ::
 
             \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
-            {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
+            {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
         if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`) or the magnetic field
+        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. 
+
+        If we write Maxwell's equations in terms of 
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
         .. math ::
 
-            \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
+            \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`).
+        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. 
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
@@ -36,7 +38,11 @@ class BaseFDEMProblem(BaseEMProblem):
 
     def fields(self, m=None):
         """
-            Solve the forward problem for the fields.
+        Solve the forward problem for the fields.
+        
+        :param numpy.array m: inversion model (nP,)
+        :rtype numpy.array:
+        :return F: forward solution
         """
 
         self.curModel = m
@@ -53,13 +59,19 @@ class BaseFDEMProblem(BaseEMProblem):
             Ainv.clean()
         return F
 
-    def Jvec(self, m, v, f=None):
+    def Jvec(self, m, v, u=None):
         """
-            Sensitivity times a vector
+        Sensitivity times a vector.
+
+        :param numpy.array m: inversion model (nP,)
+        :param numpy.array v: vector which we take sensitivity product with (nP,)
+        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :rtype numpy.array:
+        :return: Jv (ndata,) 
         """
 
-        if f is None:
-           f = self.fields(m)
+        if u is None:
+           u = self.fields(m)
 
         self.curModel = m
 
@@ -71,34 +83,41 @@ class BaseFDEMProblem(BaseEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = f[src, ftype]
+                u_src = u[src, ftype]
                 dA_dm = self.getADeriv_m(freq, u_src, v)
                 dRHS_dm = self.getRHSDeriv_m(freq, src, v) 
                 du_dm = Ainv * ( - dA_dm + dRHS_dm )
                 
                 for rx in src.rxList:
-                    df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
+                    df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
                     df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
 
-                    df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
+                    df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
                     df_dm = df_dmFun(src, v, adjoint=False)
 
+
                     Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
 
-                    P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
+                    P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
 
                     Jv[src, rx] = P(Df_Dm)
 
             Ainv.clean()
         return Utils.mkvc(Jv)
 
-    def Jtvec(self, m, v, f=None):
+    def Jtvec(self, m, v, u=None):
         """
-            Sensitivity transpose times a vector
+        Sensitivity transpose times a vector
+
+        :param numpy.array m: inversion model (nP,)
+        :param numpy.array v: vector which we take adjoint product with (nP,)
+        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :rtype numpy.array:
+        :return: Jv (ndata,) 
         """
 
-        if f is None:
-            f = self.fields(m)
+        if u is None:
+            u = self.fields(m)
 
         self.curModel = m
 
@@ -114,12 +133,12 @@ class BaseFDEMProblem(BaseEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = f[src, ftype]
+                u_src = u[src, ftype]
 
                 for rx in src.rxList:
-                    PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
 
-                    df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
+                    df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
                     df_duT = df_duTFun(src, PTv, adjoint=True)
                     
                     ATinvdf_duT = ATinv * df_duT
@@ -128,11 +147,12 @@ class BaseFDEMProblem(BaseEMProblem):
                     dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
                     du_dmT = -dA_dmT + dRHS_dmT
 
-                    df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
+                    df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
                     dfT_dm = df_dmFun(src, PTv, adjoint=True)
 
                     du_dmT += dfT_dm
 
+                    # TODO: this should be taken care of by the reciever
                     real_or_imag = rx.projComp
                     if real_or_imag is 'real':
                         Jtv +=   np.array(du_dmT,dtype=complex).real
@@ -142,15 +162,16 @@ class BaseFDEMProblem(BaseEMProblem):
                         raise Exception('Must be real or imag')
             
             ATinv.clean()
-        return Jtv
+
+        return Utils.mkvc(Jtv)
 
     def getSourceTerm(self, freq):
         """
-            Evaluates the sources for a given frequency and puts them in matrix form
+        Evaluates the sources for a given frequency and puts them in matrix form
 
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE or nF, nSrc)
-            :return: S_m, S_e
+        :param float freq: Frequency
+        :rtype: (numpy.ndarray, numpy.ndarray)  
+        :return: S_m, S_e (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._eqLocs is 'FE':
@@ -174,20 +195,22 @@ class BaseFDEMProblem(BaseEMProblem):
 
 class Problem_e(BaseFDEMProblem):
     """
-        By eliminating the magnetic flux density using
-
-            .. math ::
-
-                \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right)
-
-
-        we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+    By eliminating the magnetic flux density using
 
         .. math ::
 
-            \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right)
 
-        which we solve for \\\(\\\mathbf{e}\\\).
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+
+        \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e}
+
+    which we solve for :math:`\mathbf{e}`.
+
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _fieldType = 'e'
@@ -199,13 +222,16 @@ class Problem_e(BaseFDEMProblem):
 
     def getA(self, freq):
         """
-            .. math ::
-                \mathbf{A} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
+        System matrix
+        
+        .. math ::
+            \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
 
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
+        :param float freq: Frequency
+        :rtype: scipy.sparse.csr_matrix
+        :return: A
         """
+
         MfMui = self.MfMui
         MeSigma = self.MeSigma
         C = self.mesh.edgeCurl
@@ -214,6 +240,20 @@ class Problem_e(BaseFDEMProblem):
 
 
     def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        .. math ::
+            \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
+
+        :param float freq: frequency 
+        :param numpy.ndarray u: solution vector (nE,) 
+        :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) 
+        """
+
         dsig_dm = self.curModel.sigmaDeriv
         dMe_dsig = self.MeSigmaDeriv(u)
 
@@ -224,26 +264,37 @@ class Problem_e(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-            .. math ::
-                \mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
+        Right hand side for the system 
 
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nSrc)
-            :return: RHS
+        .. math ::
+            \mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
+
+        :param float freq: Frequency
+        :rtype: numpy.ndarray 
+        :return: RHS (nE, nSrc)
         """
 
         S_m, S_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
 
-        RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e
-
-        return RHS
+        return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
 
     def getRHSDeriv_m(self, freq, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model 
+
+        :param float freq: frequency
+        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of rhs deriv with a vector
+        """
+
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             dRHS = MfMui * (C * v)
@@ -255,20 +306,22 @@ class Problem_e(BaseFDEMProblem):
 
 class Problem_b(BaseFDEMProblem):
     """
-        We eliminate \\\(\\\mathbf{e}\\\) using
+    We eliminate :math:`\mathbf{e}` using
 
-        .. math ::
+    .. math ::
 
-             \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right)
+         \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right)
 
-        and solve for \\\(\\\mathbf{b}\\\) using:
+    and solve for :math:`\mathbf{b}` using:
 
-        .. math ::
+    .. math ::
 
-            \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}  + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
+        \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}  + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e}
 
-        .. note ::
-            The inverse problem will not work with full anisotropy
+    .. note ::
+        The inverse problem will not work with full anisotropy
+
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _fieldType = 'b'
@@ -280,12 +333,14 @@ class Problem_b(BaseFDEMProblem):
 
     def getA(self, freq):
         """
-            .. math ::
-                \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}  + i \omega
+        System matrix
 
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
+        .. math ::
+            \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}  + i \omega
+
+        :param float freq: Frequency
+        :rtype: scipy.sparse.csr_matrix
+        :return: A
         """
 
         MfMui = self.MfMui
@@ -301,6 +356,20 @@ class Problem_b(BaseFDEMProblem):
 
     def getADeriv_m(self, freq, u, v, adjoint=False):
 
+        """
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        .. math ::
+            \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
+
+        :param float freq: frequency 
+        :param numpy.ndarray u: solution vector (nF,) 
+        :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        """
+
         MfMui = self.MfMui
         C = self.mesh.edgeCurl
         MeSigmaIDeriv = self.MeSigmaIDeriv
@@ -320,12 +389,14 @@ class Problem_b(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-            .. math ::
-                \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
+        Right hand side for the system 
 
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nSrc)
-            :return: RHS
+        .. math ::
+            \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
+
+        :param float freq: Frequency
+        :rtype: numpy.ndarray 
+        :return: RHS (nE, nSrc)
         """
 
         S_m, S_e = self.getSourceTerm(freq)
@@ -341,6 +412,17 @@ class Problem_b(BaseFDEMProblem):
         return RHS
 
     def getRHSDeriv_m(self, freq, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+
+        :param float freq: frequency
+        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of rhs deriv with a vector
+        """
+
         C = self.mesh.edgeCurl
         S_m, S_e = src.eval(self)
         MfMui = self.MfMui
@@ -349,7 +431,7 @@ class Problem_b(BaseFDEMProblem):
             v = self.MfMui * v
 
         MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if not adjoint:
             RHSderiv = C * (MeSigmaIDeriv * v)
@@ -372,21 +454,22 @@ class Problem_b(BaseFDEMProblem):
 
 class Problem_j(BaseFDEMProblem):
     """
-        We eliminate \\\(\\\mathbf{h}\\\) using
+    We eliminate \\\(\\\mathbf{h}\\\) using
 
-        .. math ::
+    .. math ::
 
-            \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
+        \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
 
-        and solve for \\\(\\\mathbf{j}\\\) using
+    and solve for \\\(\\\mathbf{j}\\\) using
 
-        .. math ::
+    .. math ::
 
-            \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e}
+        \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e}
 
-        .. note::
-            This implementation does not yet work with full anisotropy!!
+    .. note::
+        This implementation does not yet work with full anisotropy!!
 
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _fieldType = 'j'
@@ -398,12 +481,14 @@ class Problem_j(BaseFDEMProblem):
 
     def getA(self, freq):
         """
-            .. math ::
-                    \\mathbf{A} = \\mathbf{C}  \\mathbf{M^e_{mu^{-1}}} \\mathbf{C}^T \\mathbf{M^f_{\\sigma^{-1}}}  + i\\omega
+        System matrix
 
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
+        .. math ::
+                \\mathbf{A} = \\mathbf{C}  \\mathbf{M^e_{\\mu^{-1}}} \\mathbf{C}^{\\top} \\mathbf{M^f_{\\sigma^{-1}}}  + i\\omega
+
+        :param float freq: Frequency
+        :rtype: scipy.sparse.csr_matrix
+        :return: A
         """
 
         MeMuI = self.MeMuI
@@ -420,12 +505,20 @@ class Problem_j(BaseFDEMProblem):
 
     def getADeriv_m(self, freq, u, v, adjoint=False):
         """
-            In this case, we assume that electrical conductivity, \\\(\\\sigma\\\) is the physical property of interest (i.e. \\\(\\\sigma\\\) = model.transform). Then we want
+        Product of the derivative of our system matrix with respect to the model and a vector
 
-            .. math ::
+        In this case, we assume that electrical conductivity, :math:`\sigma` is the physical property of interest (i.e. :math:`\sigma` = model.transform). Then we want
 
-                \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C}  \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}}
-                &= \\mathbf{C}  \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}}
+        .. math ::
+
+            \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C}  \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
+
+        :param float freq: frequency 
+        :param numpy.ndarray u: solution vector (nF,) 
+        :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
         """
 
         MeMuI = self.MeMuI
@@ -445,12 +538,15 @@ class Problem_j(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-            .. math ::
+        Right hand side for the system 
 
-                \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nSrc)
-            :return: RHS
+        .. math ::
+
+            \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
+
+        :param float freq: Frequency
+        :rtype: numpy.ndarray (nE, nSrc)
+        :return: RHS
         """
 
         S_m, S_e = self.getSourceTerm(freq)
@@ -465,9 +561,20 @@ class Problem_j(BaseFDEMProblem):
         return RHS
 
     def getRHSDeriv_m(self, freq, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model 
+
+        :param float freq: frequency
+        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of rhs deriv with a vector
+        """
+
         C = self.mesh.edgeCurl
         MeMuI = self.MeMuI
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             if self._makeASymmetric:
@@ -488,18 +595,19 @@ class Problem_j(BaseFDEMProblem):
 
 class Problem_h(BaseFDEMProblem):
     """
-        We eliminate \\\(\\\mathbf{j}\\\) using
+    We eliminate \\\(\\\mathbf{j}\\\) using
 
-        .. math ::
+    .. math ::
 
-            \mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e}
+        \mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e}
 
-        and solve for \\\(\\\mathbf{h}\\\) using
+    and solve for \\\(\\\mathbf{h}\\\) using
 
-        .. math ::
+    .. math ::
 
-            \\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
+        \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
+    :param SimPEG.Mesh mesh: mesh
     """
 
     _fieldType = 'h'
@@ -511,13 +619,14 @@ class Problem_h(BaseFDEMProblem):
 
     def getA(self, freq):
         """
-            .. math ::
+        System matrix
 
-                \mathbf{A} = \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
+        .. math::
+            \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
 
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
+        :param float freq: Frequency
+        :rtype: scipy.sparse.csr_matrix
+        :return: A
         """
 
         MeMu = self.MeMu
@@ -527,6 +636,19 @@ class Problem_h(BaseFDEMProblem):
         return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
 
     def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        .. math::
+            \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
+
+        :param float freq: frequency 
+        :param numpy.ndarray u: solution vector (nE,) 
+        :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        """
 
         MeMu = self.MeMu
         C = self.mesh.edgeCurl
@@ -538,24 +660,35 @@ class Problem_h(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-            .. math ::
+        Right hand side for the system 
 
-                \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e}
+        .. math ::
 
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, nSrc)
-            :return: RHS
+            \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
+
+        :param float freq: Frequency
+        :rtype: numpy.ndarray 
+        :return: RHS (nE, nSrc)
         """
 
         S_m, S_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
 
-        RHS = S_m + C.T * ( MfRho * S_e )
-
-        return RHS
+        return S_m + C.T * ( MfRho * S_e )
 
     def getRHSDeriv_m(self, freq, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model 
+
+        :param float freq: frequency
+        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of rhs deriv with a vector
+        """
+
         _, S_e = src.eval(self)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
@@ -566,7 +699,7 @@ class Problem_h(BaseFDEMProblem):
         elif adjoint:
             RHSDeriv = MfRhoDeriv.T * (C * v)
 
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -7,11 +7,39 @@ from SimPEG.Utils import Zero, Identity
 
 
 class Fields(SimPEG.Problem.Fields):
-    """Fancy Field Storage for a FDEM survey."""
+    """
+    
+    Fancy Field Storage for a FDEM survey. Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        b = f[srcList,'b']
+
+    If accessing all sources for a given field, use the :code:`:`
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[:,'e']
+        b = f[:,'b']
+
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
     knownFields = {}
     dtype = complex
 
 class Fields_e(Fields):
+    """
+    Fields object for Problem_e. 
+
+    :param Mesh mesh: mesh
+    :param Survey survey: survey 
+    """
+
     knownFields = {'eSolution':'E'}
     aliasFields = {
                     'e' : ['eSolution','E','_e'],
@@ -30,6 +58,15 @@ class Fields_e(Fields):
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
 
     def _ePrimary(self, eSolution, srcList):
+        """
+        Primary electric field from source
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary electric field as defined by the sources
+        """
+
         ePrimary = np.zeros_like(eSolution)
         for i, src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
@@ -37,19 +74,67 @@ class Fields_e(Fields):
         return ePrimary
 
     def _eSecondary(self, eSolution, srcList):
+        """
+        Secondary electric field is the thing we solved for
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary electric field
+        """
+
         return eSolution
 
     def _e(self, eSolution, srcList):
+        """
+        Total electric field is sum of primary and secondary 
+        
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total electric field
+        """
+
         return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
 
     def _eDeriv_u(self, src, v, adjoint = False):
+        """
+        Derivative of the total electric field with respect to the thing we 
+        solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the electric field with respect to the field we solved for with a vector
+        """
+
         return Identity()*v
 
     def _eDeriv_m(self, src, v, adjoint = False):
+        """
+        Derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the electric field derivative with respect to the inversion model with a vector
+        """
+
         # assuming primary does not depend on the model
         return Zero()
 
     def _bPrimary(self, eSolution, srcList):
+        """
+        Primary magnetic flux density from source
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density as defined by the sources
+        """
+
         bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
@@ -57,6 +142,15 @@ class Fields_e(Fields):
         return bPrimary
 
     def _bSecondary(self, eSolution, srcList):
+        """
+        Secondary magnetic flux density from eSolution
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary magnetic flux density
+        """
+
         C = self._edgeCurl
         b = (C * eSolution)
         for i, src in enumerate(srcList):
@@ -66,29 +160,84 @@ class Fields_e(Fields):
         return b
 
     def _bSecondaryDeriv_u(self, src, v, adjoint = False):
+        """
+        Derivative of the secondary magnetic flux density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector
+        """
+
         C = self._edgeCurl
         if adjoint:
             return - 1./(1j*omega(src.freq)) * (C.T * v)
         return - 1./(1j*omega(src.freq)) * (C * v)
 
     def _bSecondaryDeriv_m(self, src, v, adjoint = False):
-        S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
-        S_mDeriv = S_mDeriv(v)
+        """
+        Derivative of the secondary magnetic flux density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector
+        """
+
+        S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
         return 1./(1j * omega(src.freq)) * S_mDeriv
 
     def _b(self, eSolution, srcList):
+        """
+        Total magnetic flux density is sum of primary and secondary 
+        
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic flux density 
+        """
+
         return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
 
     def _bDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic flux density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
+        """
+
         # Primary does not depend on u
         return self._bSecondaryDeriv_u(src, v, adjoint)
 
     def _bDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic flux density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
+        """
+
         # Assuming the primary does not depend on the model
         return self._bSecondaryDeriv_m(src, v, adjoint)
 
 
 class Fields_b(Fields):
+    """
+    Fields object for Problem_b. 
+
+    :param Mesh mesh: mesh
+    :param Survey survey: survey 
+    """
+
     knownFields = {'bSolution':'F'}
     aliasFields = {
                     'b' : ['bSolution','F','_b'],
@@ -111,6 +260,15 @@ class Fields_b(Fields):
         self._Me = self.survey.prob.Me
 
     def _bPrimary(self, bSolution, srcList):
+        """
+        Primary magnetic flux density from source
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary electric field as defined by the sources
+        """
+
         bPrimary = np.zeros_like(bSolution)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
@@ -118,19 +276,66 @@ class Fields_b(Fields):
         return bPrimary
 
     def _bSecondary(self, bSolution, srcList):
+        """
+        Secondary magnetic flux density is the thing we solved for
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary magnetic flux density
+        """
+
         return bSolution
 
     def _b(self, bSolution, srcList):
+        """
+        Total magnetic flux density is sum of primary and secondary 
+        
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic flux density 
+        """
+
         return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
 
     def _bDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic flux density with respect to the thing we 
+        solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
+        """
         return Identity()*v
 
     def _bDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
+        """
+
         # assuming primary does not depend on the model
         return Zero()
 
     def _ePrimary(self, bSolution, srcList):
+        """
+        Primary electric field from source
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary electric field as defined by the sources
+        """
+
         ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
         for i,src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
@@ -138,6 +343,15 @@ class Fields_b(Fields):
         return ePrimary
 
     def _eSecondary(self, bSolution, srcList):
+        """
+        Secondary electric field from bSolution
+
+        :param numpy.ndarray bSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary electric field
+        """
+
         e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
@@ -145,12 +359,32 @@ class Fields_b(Fields):
         return e
 
     def _eSecondaryDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the secondary electric field with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary electric field with respect to the field we solved for with a vector
+        """
+
         if not adjoint:
             return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) )
         else:
             return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
 
     def _eSecondaryDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the secondary electric field with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary electric field with respect to the model with a vector
+        """
+
         bSolution = self[[src],'bSolution']
         _,S_e = src.eval(self.prob)
         Me = self._Me
@@ -166,25 +400,60 @@ class Fields_b(Fields):
         elif adjoint:
             de_dm = self._MeSigmaIDeriv(w).T * v
 
-        _, S_eDeriv = src.evalDeriv(self.prob, adjoint)
-        Se_Deriv = S_eDeriv(v)
+        _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
 
-        de_dm = de_dm - self._MeSigmaI * Se_Deriv
+        de_dm = de_dm - self._MeSigmaI * S_eDeriv
 
         return de_dm
 
     def _e(self, bSolution, srcList):
+        """
+        Total electric field is sum of primary and secondary 
+        
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total electric field
+        """
+
         return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList)
 
     def _eDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total electric field with respect to the thing we solved for 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the electric field with respect to the field we solved for with a vector
+        """
+
         return self._eSecondaryDeriv_u(src, v, adjoint)
 
     def _eDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total electric field density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the electric field derivative with respect to the inversion model with a vector
+        """
+
         # assuming primary doesn't depend on model
         return self._eSecondaryDeriv_m(src, v, adjoint)
 
 
 class Fields_j(Fields):
+    """
+    Fields object for Problem_j. 
+
+    :param Mesh mesh: mesh
+    :param Survey survey: survey 
+    """
+
     knownFields = {'jSolution':'F'}
     aliasFields = {
                     'j' : ['jSolution','F','_j'],
@@ -207,6 +476,15 @@ class Fields_j(Fields):
         self._Me = self.survey.prob.Me
 
     def _jPrimary(self, jSolution, srcList):
+        """
+        Primary current density from source
+
+        :param numpy.ndarray jSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary current density as defined by the sources
+        """
+
         jPrimary = np.zeros_like(jSolution,dtype = complex)
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
@@ -214,19 +492,66 @@ class Fields_j(Fields):
         return jPrimary
 
     def _jSecondary(self, jSolution, srcList):
+        """
+        Secondary current density is the thing we solved for
+
+        :param numpy.ndarray jSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary current density
+        """
+
         return jSolution
 
     def _j(self, jSolution, srcList):
+        """
+        Total current density is sum of primary and secondary 
+        
+        :param numpy.ndarray jSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total current density 
+        """
+
         return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
 
     def _jDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total current density with respect to the thing we 
+        solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the current density with respect to the field we solved for with a vector
+        """
+
         return Identity()*v
 
     def _jDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the current density derivative with respect to the inversion model with a vector
+        """
         # assuming primary does not depend on the model
         return Zero()
 
     def _hPrimary(self, jSolution, srcList):
+        """
+        Primary magnetic field from source
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary magnetic field as defined by the sources
+        """
+
         hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
             hp = src.hPrimary(self.prob)
@@ -234,6 +559,15 @@ class Fields_j(Fields):
         return hPrimary
 
     def _hSecondary(self, jSolution, srcList):
+        """
+        Secondary magnetic field from bSolution
+
+        :param numpy.ndarray jSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary magnetic field
+        """
+
         h =  self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) )
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
@@ -242,12 +576,32 @@ class Fields_j(Fields):
         return h
 
     def _hSecondaryDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the secondary magnetic field with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector
+        """
+
         if not adjoint:
             return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) )
         elif adjoint:
             return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v))
 
     def _hSecondaryDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the secondary magnetic field with respect to the inversion model 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary magnetic field with respect to the model with a vector
+        """
+
         jSolution = self[[src],'jSolution']
         MeMuI = self._MeMuI
         C = self._edgeCurl
@@ -260,7 +614,7 @@ class Fields_j(Fields):
         elif adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
 
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
+        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
         if not adjoint:
             S_mDeriv = S_mDeriv(v)
@@ -272,17 +626,53 @@ class Fields_j(Fields):
 
 
     def _h(self, jSolution, srcList):
+        """
+        Total magnetic field is sum of primary and secondary 
+        
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic field
+        """
+
         return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
 
     def _hDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic field with respect to the thing we solved for 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
+        """
+
         return self._hSecondaryDeriv_u(src, v, adjoint)
 
     def _hDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic field density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the magnetic field derivative with respect to the inversion model with a vector
+        """
+
         # assuming the primary doesn't depend on the model
         return self._hSecondaryDeriv_m(src, v, adjoint)
 
 
 class Fields_h(Fields):
+    """
+    Fields object for Problem_h. 
+
+    :param Mesh mesh: mesh
+    :param Survey survey: survey 
+    """
+
     knownFields = {'hSolution':'E'}
     aliasFields = {
                     'h' : ['hSolution','E','_h'],
@@ -303,6 +693,15 @@ class Fields_h(Fields):
         self._MfRho = self.survey.prob.MfRho
 
     def _hPrimary(self, hSolution, srcList):
+        """
+        Primary magnetic field from source
+
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary magnetic field as defined by the sources
+        """
+
         hPrimary = np.zeros_like(hSolution,dtype = complex)
         for i, src in enumerate(srcList):
             hp = src.hPrimary(self.prob)
@@ -310,19 +709,67 @@ class Fields_h(Fields):
         return hPrimary
 
     def _hSecondary(self, hSolution, srcList):
+        """
+        Secondary magnetic field is the thing we solved for
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary magnetic field
+        """
+
         return hSolution
 
     def _h(self, hSolution, srcList):
+        """
+        Total magnetic field is sum of primary and secondary 
+        
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic field
+        """
+
         return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
 
     def _hDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic field with respect to the thing we 
+        solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
+        """
+
         return Identity()*v
 
     def _hDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the magnetic field derivative with respect to the inversion model with a vector
+        """
+
         # assuming primary does not depend on the model
         return Zero()
 
     def _jPrimary(self, hSolution, srcList):
+        """
+        Primary current density from source
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: primary current density as defined by the sources
+        """
+
         jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex)
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
@@ -330,6 +777,15 @@ class Fields_h(Fields):
         return jPrimary
 
     def _jSecondary(self, hSolution, srcList):
+        """
+        Secondary current density from eSolution
+
+        :param numpy.ndarray hSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: secondary current density
+        """
+
         j = self._edgeCurl*hSolution
         for i, src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
@@ -337,22 +793,69 @@ class Fields_h(Fields):
         return j
 
     def _jSecondaryDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the secondary current density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector
+        """
+
         if not adjoint:
             return self._edgeCurl*v
         elif adjoint:
             return self._edgeCurl.T*v
 
     def _jSecondaryDeriv_m(self, src, v, adjoint=False):
-        _,S_eDeriv = src.evalDeriv(self.prob, adjoint)
-        S_eDeriv = S_eDeriv(v)
+        """
+        Derivative of the secondary current density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the secondary current density derivative with respect to the inversion model with a vector
+        """
+
+        _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
         return -S_eDeriv
 
     def _j(self, hSolution, srcList):
+        """
+        Total current density is sum of primary and secondary 
+        
+        :param numpy.ndarray eSolution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total current density 
+        """
+
         return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
 
     def _jDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the total current density with respect to the thing we solved for
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of the derivative of the current density with respect to the field we solved for with a vector
+        """
         return self._jSecondaryDeriv_u(src,v,adjoint)
 
     def _jDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the total current density with respect to the inversion model. 
+        
+        :param SimPEG.EM.FDEM.Src src: source
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: SimPEG.Utils.Zero
+        :return: product of the current density with respect to the inversion model with a vector
+        """
+
         # assuming the primary does not depend on the model
         return self._jSecondaryDeriv_m(src,v,adjoint)

SimPEG/EM/FDEM/SrcFDEM.py

@@ -1,55 +1,141 @@
 from SimPEG import Survey, Problem, Utils, np, sp
 from scipy.constants import mu_0
 from SimPEG.EM.Utils import *
-from SimPEG.Utils import Zero
-# from SurveyFDEM import Rx
- 
+from SimPEG.Utils import Zero 
 
 class BaseSrc(Survey.BaseSrc):
+    """
+    Base source class for FDEM Survey
+    """
+
     freq = None
-    # rxPair = Rx
+    # rxPair = RxFDEM
     integrate = True
 
     def eval(self, prob):
+        """
+        Evaluate the source terms. 
+        - :math:`S_m` : magnetic source term
+        - :math:`S_e` : electric source term
+
+        :param Problem prob: FDEM Problem
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: tuple with magnetic source term and electric source term
+        """
         S_m = self.S_m(prob)
         S_e = self.S_e(prob)
         return S_m, S_e
 
-    def evalDeriv(self, prob, v, adjoint=False):
-        return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
+    def evalDeriv(self, prob, v=None, adjoint=False):
+        """
+        Derivatives of the source terms with respect to the inversion model
+        - :code:`S_mDeriv` : derivative of the magnetic source term
+        - :code:`S_eDeriv` : derivative of the electric source term
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: tuple with magnetic source term and electric source term derivatives times a vector 
+        """
+        if v is not None: 
+            return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint)
+        else:
+            return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
 
     def bPrimary(self, prob):
+        """
+        Primary magnetic flux density
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density 
+        """
         return Zero()
 
     def hPrimary(self, prob):
+        """
+        Primary magnetic field
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         return Zero()
 
     def ePrimary(self, prob):
+        """
+        Primary electric field
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: primary electric field 
+        """
         return Zero()
 
     def jPrimary(self, prob):
+        """
+        Primary current density
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: primary current density 
+        """
         return Zero()
 
     def S_m(self, prob):
+        """
+        Magnetic source term 
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: magnetic source term on mesh
+        """
         return Zero()
 
     def S_e(self, prob):
+        """
+        Electric source term 
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: electric source term on mesh
+        """
         return Zero()
 
     def S_mDeriv(self, prob, v, adjoint = False):
+        """
+        Derivative of magnetic source term with respect to the inversion model
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of magnetic source term derivative with a vector
+        """
+
         return Zero()
 
     def S_eDeriv(self, prob, v, adjoint = False):
+        """
+        Derivative of electric source term with respect to the inversion model
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector to take product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: product of electric source term derivative with a vector
+        """
         return Zero()
 
 
 class RawVec_e(BaseSrc):
     """
-        RawVec electric source. It is defined by the user provided vector S_e
+    RawVec electric source. It is defined by the user provided vector S_e
 
-        :param numpy.array S_e: electric source term
-        :param float freq: frequency
-        :param rxList: receiver list
+    :param list rxList: receiver list
+    :param float freq: frequency
+    :param numpy.array S_e: electric source term
     """
 
     def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
@@ -58,16 +144,17 @@ class RawVec_e(BaseSrc):
         BaseSrc.__init__(self, rxList)
 
     def S_e(self, prob):
+
         return self._S_e
 
 
 class RawVec_m(BaseSrc):
     """
-        RawVec magnetic source. It is defined by the user provided vector S_m
+    RawVec magnetic source. It is defined by the user provided vector S_m
 
-        :param numpy.array S_m: magnetic source term
-        :param float freq: frequency
-        :param rxList: receiver list
+    :param float freq: frequency
+    :param rxList: receiver list
+    :param numpy.array S_m: magnetic source term
     """
 
     def __init__(self, rxList, freq, S_m, integrate = True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
@@ -78,17 +165,24 @@ class RawVec_m(BaseSrc):
         BaseSrc.__init__(self, rxList)
 
     def S_m(self, prob):
+        """
+        Magnetic source term 
+
+        :param Problem prob: FDEM Problem
+        :rtype: numpy.ndarray
+        :return: magnetic source term on mesh
+        """
         return self._S_m
 
 
 class RawVec(BaseSrc):
     """
-        RawVec source. It is defined by the user provided vectors S_m, S_e
+    RawVec source. It is defined by the user provided vectors S_m, S_e
 
-        :param numpy.array S_m: magnetic source term
-        :param numpy.array S_e: electric source term
-        :param float freq: frequency
-        :param rxList: receiver list
+    :param rxList: receiver list
+    :param float freq: frequency
+    :param numpy.array S_m: magnetic source term
+    :param numpy.array S_e: electric source term
     """
     def __init__(self, rxList, freq, S_m, S_e, integrate = True):
         self._S_m = np.array(S_m,dtype=complex)
@@ -109,6 +203,51 @@ class RawVec(BaseSrc):
 
 
 class MagDipole(BaseSrc):
+    """  
+    Point magnetic dipole source calculated by taking the curl of a magnetic
+    vector potential. By taking the discrete curl, we ensure that the magnetic
+    flux density is divergence free (no magnetic monopoles!). 
+
+    This approach uses a primary-secondary in frequency. Here we show the
+    derivation for E-B formulation noting that similar steps are followed for
+    the H-J formulation.
+
+    .. math:: 
+        \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
+            {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
+
+    We split up the fields and :math:`\mu^{-1}` into primary (:math:`\mathbf{P}`) and secondary (:math:`\mathbf{S}`) components
+
+    - :math:`\mathbf{e} = \mathbf{e^P} + \mathbf{e^S}`
+    - :math:`\mathbf{b} = \mathbf{b^P} + \mathbf{b^S}`
+    - :math:`\\boldsymbol{\mu}^{\mathbf{-1}} = \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{P}} + \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{S}}`
+
+    and define a zero-frequency primary problem, noting that the source is
+    generated by a divergence free electric current
+
+    .. math:: 
+        \mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\
+            {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}}
+
+    Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is 
+
+    .. math:: 
+        \mathbf{e^P} =  0 \\\\
+            {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}}
+
+    Our secondary problem is then 
+
+    .. math::
+        \mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\
+            {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b^S} - \mathbf{M_{\sigma}^e} \mathbf{e^S} = -\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^S} \mathbf{b^P}}
+
+    :param list rxList: receiver list
+    :param float freq: frequency
+    :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
+    :param string orientation: 'X', 'Y', 'Z'
+    :param float moment: magnetic dipole moment 
+    :param float mu: background magnetic permeability 
+    """
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
     def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
@@ -121,6 +260,13 @@ class MagDipole(BaseSrc):
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
+        """
+        The primary magnetic flux density from a magnetic vector potential
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -152,14 +298,37 @@ class MagDipole(BaseSrc):
         return C*a
 
     def hPrimary(self, prob):
+        """
+        The primary magnetic field from a magnetic vector potential
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         b = self.bPrimary(prob)
         return h_from_b(prob,b)
 
     def S_m(self, prob):
+        """
+        The magnetic source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
+
         b_p = self.bPrimary(prob)
         return -1j*omega(self.freq)*b_p
 
     def S_e(self, prob):
+        """
+        The electric source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
+
         if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
             return Zero()
         else:
@@ -179,6 +348,21 @@ class MagDipole(BaseSrc):
 
 class MagDipole_Bfield(BaseSrc):
 
+    """
+    Point magnetic dipole source calculated with the analytic solution for the
+    fields from a magnetic dipole. No discrete curl is taken, so the magnetic
+    flux density may not be strictly divergence free.
+
+    This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. 
+
+    :param list rxList: receiver list
+    :param float freq: frequency
+    :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
+    :param string orientation: 'X', 'Y', 'Z'
+    :param float moment: magnetic dipole moment 
+    :param float mu: background magnetic permeability 
+    """
+
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
     #TODO: neither does moment
     def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
@@ -190,6 +374,14 @@ class MagDipole_Bfield(BaseSrc):
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
+        """
+        The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
+
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -221,14 +413,35 @@ class MagDipole_Bfield(BaseSrc):
         return b
 
     def hPrimary(self, prob):
+        """
+        The primary magnetic field from a magnetic vector potential
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         b = self.bPrimary(prob)
         return h_from_b(prob, b)
 
     def S_m(self, prob):
+        """
+        The magnetic source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         b = self.bPrimary(prob)
         return -1j*omega(self.freq)*b
 
     def S_e(self, prob):
+        """
+        The electric source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
             return Zero()
         else:
@@ -247,6 +460,20 @@ class MagDipole_Bfield(BaseSrc):
 
 
 class CircularLoop(BaseSrc):
+    """
+    Circular loop magnetic source calculated by taking the curl of a magnetic
+    vector potential. By taking the discrete curl, we ensure that the magnetic
+    flux density is divergence free (no magnetic monopoles!). 
+
+    This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. 
+
+    :param list rxList: receiver list
+    :param float freq: frequency
+    :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
+    :param string orientation: 'X', 'Y', 'Z'
+    :param float moment: magnetic dipole moment 
+    :param float mu: background magnetic permeability 
+    """
 
     #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
     def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
@@ -259,6 +486,13 @@ class CircularLoop(BaseSrc):
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
+        """
+        The primary magnetic flux density from a magnetic vector potential
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         eqLocs = prob._eqLocs
 
         if eqLocs is 'FE':
@@ -277,10 +511,10 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', self.radius, mu=self.mu)
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticLoopVectorPotential
+            srcfct = MagneticDipoleVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -289,14 +523,35 @@ class CircularLoop(BaseSrc):
         return C*a
 
     def hPrimary(self, prob):
+        """
+        The primary magnetic field from a magnetic vector potential
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         b = self.bPrimary(prob)
         return 1./self.mu*b
 
     def S_m(self, prob):
+        """
+        The magnetic source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         b = self.bPrimary(prob)
         return -1j*omega(self.freq)*b
 
     def S_e(self, prob):
+        """
+        The electric source term
+
+        :param Problem prob: FDEM problem
+        :rtype: numpy.ndarray
+        :return: primary magnetic field 
+        """
         if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
             return Zero()
         else:

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -10,6 +10,12 @@ import SrcFDEM as Src
 ####################################################
 
 class Rx(SimPEG.Survey.BaseRx):
+    """
+    Frequency domain receivers
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string rxType: reciever type from knownRxTypes
+    """
 
     knownRxTypes = {
                     'exr':['e', 'Ex', 'real'],
@@ -26,13 +32,6 @@ class Rx(SimPEG.Survey.BaseRx):
                     'byi':['b', 'Fy', 'imag'],
                     'bzi':['b', 'Fz', 'imag'],
 
-                    'bxr_sec':['bSecondary', 'Fx', 'real'],
-                    'byr_sec':['bSecondary', 'Fy', 'real'],
-                    'bzr_sec':['bSecondary', 'Fz', 'real'],
-                    'bxi_sec':['bSecondary', 'Fx', 'imag'],
-                    'byi_sec':['bSecondary', 'Fy', 'imag'],
-                    'bzi_sec':['bSecondary', 'Fz', 'imag'],
-
                     'jxr':['j', 'Fx', 'real'],
                     'jyr':['j', 'Fy', 'real'],
                     'jzr':['j', 'Fz', 'real'],
@@ -68,6 +67,15 @@ class Rx(SimPEG.Survey.BaseRx):
         return self.knownRxTypes[self.rxType][2]
 
     def projectFields(self, src, mesh, u):
+        """
+        Project fields to recievers to get data.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
         P = self.getP(mesh)
         u_part_complex = u[src, self.projField]
         # get the real or imag component
@@ -76,6 +84,16 @@ class Rx(SimPEG.Survey.BaseRx):
         return P*u_part
 
     def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
+        """
+        Derivative of projected fields with respect to the inversion model times a vector.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields u: fields object
+        :param numpy.ndarray v: vector to multiply
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
         P = self.getP(mesh)
 
         if not adjoint:
@@ -102,10 +120,13 @@ class Rx(SimPEG.Survey.BaseRx):
 
 class Survey(SimPEG.Survey.BaseSurvey):
     """
-        docstring for SurveyFDEM
+    Frequency domain electromagnetic survey
+
+    :param list srcList: list of FDEM sources used in the survey
     """
 
     srcPair = Src.BaseSrc
+    rxPaair = Rx 
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
@@ -113,7 +134,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
         SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
         _freqDict = {}
-        for src in self.srcList:
+        for src in srcList:
             if src.freq not in _freqDict:
                 _freqDict[src.freq] = []
             _freqDict[src.freq] += [src]
@@ -133,6 +154,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
 
     @property
     def nSrcByFreq(self):
+        """Number of sources at each frequency"""
         if getattr(self, '_nSrcByFreq', None) is None:
             self._nSrcByFreq = {}
             for freq in self.freqs:
@@ -140,11 +162,22 @@ class Survey(SimPEG.Survey.BaseSurvey):
         return self._nSrcByFreq
 
     def getSrcByFreq(self, freq):
-        """Returns the sources associated with a specific frequency."""
+        """
+        Returns the sources associated with a specific frequency.
+        :param float freq: frequency for which we look up sources
+        :rtype: dictionary
+        :return: sources at the sepcified frequency 
+        """
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
     def projectFields(self, u):
+        """
+        Project fields to receiver locations
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: data
+        """
         data = SimPEG.Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -79,32 +79,12 @@ class SrcTDEM(Survey.BaseSrc):
 
 class SrcTDEM_VMD_MVP(SrcTDEM):
 
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+    def __init__(self,rxList,loc):
         self.loc = loc
-        self.waveformType = waveformType
         SrcTDEM.__init__(self,rxList)
 
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"        
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD') 
-            return {"b": mesh.edgeCurl*MVP}               
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
                 MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -113,12 +93,13 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
-            raise Exception('Unknown mesh for VMD')            
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+            raise Exception('Unknown mesh for VMD')
+
+        return {"b": mesh.edgeCurl*MVP}
 
 
 class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+    def __init__(self,rxList,loc,radius,waveformType):
         self.loc = loc
         self.radius = radius
         self.waveformType = waveformType

SimPEG/EM/__init__.py

@@ -1,6 +1,6 @@
-# from EM import *
 import TDEM
 import FDEM
 import Base
 import Analytics
 import Utils
+from scipy.constants import mu_0, epsilon_0

SimPEG/Examples/DC_PseudoSection_Simulation.py

@@ -1,179 +0,0 @@
-from SimPEG import *
-import simpegDCIP as DC
-import scipy.interpolate as interpolation
-import matplotlib.pyplot as plt
-import time
-import re
-
-def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp',  plotIt=True):
-    """
-    DC Forward Simulation
-    
-    Forward model conductive spheres in a half-space and plot a pseudo-section
-    
-    Created on Mon Feb 01 19:28:06 2016
-    
-    @fourndo
-    """
-    
-    # First we need to create a mesh and a model.
-    
-    # This is our mesh
-    dx    = 5.
-    
-    hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)]
-    hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)]
-    hzind = [(dx,15,-1.3),(dx, 15)]
-    
-    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
-    
-    
-    # Set background conductivity
-    model = np.ones(mesh.nC) * sig[0]
-    
-    # First anomaly
-    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
-    model[ind] = sig[1]
-    
-    # Second anomaly
-    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC)
-    model[ind] = sig[2]
-    
-    # Get index of the center
-    indy = int(mesh.nCy/2)
-
-    
-    # Plot the model for reference
-    # Define core mesh extent
-    xlim = 200
-    zlim = 125
-    
-    # Specify the survey type: "pdp" | "dpdp"
-    
-    
-    # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
-    ends = [(-175,0),(175,0)]
-    ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]]
-    
-    # Snap the endpoints to the grid. Easier to create 2D section.
-    indx = Utils.closestPoints(mesh, ends )
-    locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
-    
-    # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
-    [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
-    
-    # Define some global geometry
-    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) 
-    dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
-    dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
-    azm =  np.arctan(dl_y/dl_x)
-
-    #Set boundary conditions
-    mesh.setCellGradBC('neumann')
-    
-    # Define the differential operators needed for the DC problem
-    Div = mesh.faceDiv
-    Grad = mesh.cellGrad
-    Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
-    
-    A = Div*Msig*Grad
-    
-    # Change one corner to deal with nullspace
-    A[0,0] = 1
-    A = sp.csc_matrix(A)
-    
-    # We will solve the system iteratively, so a pre-conditioner is helpful
-    # This is simply a Jacobi preconditioner (inverse of the main diagonal)
-    dA = A.diagonal()
-    P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
-
-    # Now we can solve the system for all the transmitters
-    # We want to store the data
-    data = []
-    
-    # There is probably a more elegant way to do this, but we can just for-loop through the transmitters
-    for ii in range(len(Tx)):
-        
-        start_time = time.time() # Let's time the calculations
-        
-        #print("Transmitter %i / %i\r" % (ii+1,len(Tx)))
-        
-        # Select dipole locations for receiver
-        rxloc_M = np.asarray(Rx[ii][:,0:3])
-        rxloc_N = np.asarray(Rx[ii][:,3:])
-
-      
-        # For usual cases "dpdp" or "gradient"
-        if not re.match(stype,'pdp'): 
-            inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )   
-            
-        else: 
-            
-            # Create an "inifinity" pole
-            tx =  np.squeeze(Tx[ii][:,0:1])
-            tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
-            inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )  
-    
-    
-        # Iterative Solve
-        Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
-    
-        # We now have the potential everywhere
-        phi = mkvc(Ainvb[0])
-        
-        # Solve for phi on pole locations
-        P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
-        P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
-            
-        # Compute the potential difference
-        dtemp = (P1*phi - P2*phi)*np.pi
-        
-        data.append( dtemp )     
-        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
-        
-    print 'Transmitter {0} of {1}'.format(ii,len(Tx))
-    print 'Forward completed'
-    
-    
-    # Let's just convert the 3D format into 2D (distance along line) and plot
-    [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx)
-    
-
-    # Here is an example for the first tx-rx array
-    if plotIt:
-        fig = plt.figure()
-        ax = plt.subplot(2,1,1, aspect='equal')
-        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
-        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
-        plt.gca().set_aspect('equal', adjustable='box')
-
-        plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
-        plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
-        plt.xlim([-xlim,xlim])
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
-        
-        
-        ax = plt.subplot(2,1,2, aspect='equal')   
-        
-        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax.add_artist(circle1)
-        ax.add_artist(circle2)  
-        
-        # Add the speudo section
-        DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype)
-        
-        plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
-        plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
-        plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
-        
-        plt.show()
-        
-        return fig, ax
-        
-if __name__ == '__main__':
-    run()

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -0,0 +1,116 @@
+from SimPEG import *
+import SimPEG.EM as EM
+from SimPEG.EM import mu_0
+
+
+def run(plotIt=True):
+    """
+        EM: FDEM: 1D: Inversion
+        =======================
+
+        Here we will create and run a FDEM 1D inversion.
+
+    """
+
+    cs, ncx, ncz, npad = 5., 25, 15, 15
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+    layerz = -100.
+
+    active = mesh.vectorCCz<0.
+    layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
+    actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
+    sig_half = 2e-2
+    sig_air = 1e-8
+    sig_layer = 1e-2
+    sigma = np.ones(mesh.nCz)*sig_air
+    sigma[active] = sig_half
+    sigma[layer] = sig_layer
+    mtrue = np.log(sigma[active])
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (3, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active])
+        ax.set_ylim(-500, 0)
+        ax.set_xlim(1e-3, 1e-1)
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+
+
+    rxOffset=10. 
+    bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
+
+    freqs = np.logspace(1,3,10)
+    srcLoc = np.array([0., 0., 10.])
+
+    srcList = []
+    [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+
+    survey = EM.FDEM.Survey(srcList)
+    prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
+
+    prb.pair(survey)
+
+    std = 0.05
+    survey.makeSyntheticData(mtrue, std)
+
+    survey.std = std
+    survey.eps = np.linalg.norm(survey.dtrue)*1e-5
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (6, 6))
+        ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-')
+        ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-')
+        ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
+        ax.set_xlabel('Time (s)', fontsize = 14)
+        ax.set_ylabel('$B_z$ (T)', fontsize = 16)
+        ax.set_xlabel('Time (s)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+
+    dmisfit = DataMisfit.l2_DataMisfit(survey)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Tikhonov(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIter = 6)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
+    # Create an inversion object
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
+    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 1e-3
+    reg.alpha_x = 1.
+    prb.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+
+    mopt = inv.run(m0)
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (3, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active])
+        plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
+        ax.set_ylim(-500, 0)
+        ax.set_xlim(1e-3, 1e-1)
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+        plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best')
+        plt.show()
+
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_FDEM_SusEffects.py

@@ -0,0 +1,154 @@
+from SimPEG import *
+from SimPEG import EM
+from scipy.constants import mu_0
+
+def run(plotIt=True):
+    """
+        EM: FDEM: Effects of susceptibility
+        ===================================
+
+        When airborne freqeuncy domain EM (AFEM) survey is flown over
+        the earth including significantly susceptible bodies (magnetite-rich rocks),
+        negative data is often observed in the real part of the lowest frequency
+        (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
+        occurs due to a susceptible body when the magnetic field is applied.
+
+        To clarify what is happening in the earth when we are exciting the earth with
+        a loop source in the frequency domain we run three forward modelling:
+
+            - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
+            - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
+            - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
+
+        We plot vector magnetic fields in the earth. For secondary fields we provide
+        F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
+        figure show both real and parts.
+
+    """
+    # Generate Cylindrical mesh
+    cs, ncx, ncz, npad = 5, 25, 24, 20.
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+    sighalf = 1e-3
+    sigma = np.ones(mesh.nC)*1e-8
+    sigmahomo = sigma.copy()
+    mu = np.ones(mesh.nC)*mu_0
+    sigma[mesh.gridCC[:,-1]<0.] = sighalf
+    blkind  = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50))
+    sigma[blkind] = 1e-1
+    mu[blkind] = mu_0*1.1
+    offset = 0.
+    frequency = np.r_[10., 100., 1000.]
+    rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr')
+    rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi')
+    srcLists = []
+    nfreq = frequency.size
+    for ifreq in range(nfreq):
+        src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
+        srcLists.append(src)
+    survey = EM.FDEM.Survey(srcLists)
+    iMap = Maps.IdentityMap(nP=int(mesh.nC))
+    # Use PhysPropMap
+    maps = [('sigma', iMap), ('mu', iMap)]
+    prob = EM.FDEM.Problem_b(mesh, mapping=maps)
+    prob.Solver = MumpsSolver
+    survey.pair(prob)
+    m = np.r_[sigma, mu]
+    survey0 = EM.FDEM.Survey(srcLists)
+    prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
+    
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+            prb.Solver = SolverLU
+
+    survey0.pair(prob0)
+    m = np.r_[sigma, mu]
+    m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
+    m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
+    # Anomalous conductivity and susceptibility
+    F = prob.fields(m)
+    # Only anomalous conductivity
+    F0 = prob.fields(m0)
+    # Primary field
+    F00 = prob.fields(m00)
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        def vizfields(ifreq=0, primsec="secondary",realimag="real"):
+
+            titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"]
+            actind  = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200))
+
+            if primsec=="secondary":
+                bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
+                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
+            elif primsec=="primary":
+                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+
+            XYZ = mesh.gridCC[actind,:]
+            X = XYZ[:,0].reshape((31,43), order='F')
+            Z = XYZ[:,2].reshape((31,43), order='F')
+            bx = bCC[actind,0].reshape((31,43), order='F')
+            bz = bCC[actind,1].reshape((31,43), order='F')
+            bx0 = bCC0[actind,0].reshape((31,43), order='F')
+            bz0 = bCC0[actind,1].reshape((31,43), order='F')
+
+            bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F')
+            bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F')
+
+            absbreal = np.sqrt(bx.real**2+bz.real**2)
+            absbimag = np.sqrt(bx.imag**2+bz.imag**2)
+            absb0real = np.sqrt(bx0.real**2+bz0.real**2)
+            absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2)
+
+            absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2)
+            absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2)
+
+            fig = plt.figure(figsize=(15,5))
+            ax1 = plt.subplot(131)
+            ax2 = plt.subplot(132)
+            ax3 = plt.subplot(133)
+            typefield="real"
+            scale=20
+            if realimag=="real":
+                ax1.contourf(X, Z,np.log10(absbreal), 100)
+                ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5)
+                ax2.contourf(X, Z,np.log10(absb0real), 100)
+                ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5)
+                ax3.contourf(X, Z,np.log10(absbrealsec), 100)
+                ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5)
+            elif realimag=="imag":
+                ax1.contourf(X, Z,np.log10(absbimag), 100)
+                ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5)
+                ax2.contourf(X, Z,np.log10(absb0imag), 100)
+                ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5)
+                ax3.contourf(X, Z,np.log10(absbimagsec), 100)
+                ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5)
+
+            ax = [ax1, ax2, ax3]
+            ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18)
+            ax2.text(30, 50, primsec, color="k", fontsize=18)
+            ax1.text(30, 50, realimag, color="k", fontsize=18)
+            for i, axtemp in enumerate(ax):
+                axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3)
+
+                axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3)
+                axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3)
+                axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3)
+                axtemp.set_ylim(-200, 100.)
+                axtemp.set_xlim(10, 100.)
+                axtemp.set_title(titles[i])
+            plt.show()
+            return fig, ax
+        fig1, ax1 = vizfields(1, primsec="primary", realimag="real")
+        fig2, ax2 = vizfields(1, primsec="secondary", realimag="real")
+        fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag")
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -1,6 +1,6 @@
 from SimPEG import *
 import SimPEG.EM as EM
-from scipy.constants import mu_0
+from SimPEG.EM import mu_0
 
 
 def run(plotIt=True):
@@ -50,20 +50,18 @@ def run(plotIt=True):
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
     prb.pair(survey)
-    dtrue = survey.dpred(mtrue)
 
-
-    survey.dtrue = dtrue
+    # create observed data
     std = 0.05
-    noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
-    survey.dobs = survey.dtrue+noise
-    survey.std = survey.dobs*0 + std
-    survey.Wd = 1/(abs(survey.dobs)*std)
+    
+    survey.dobs = survey.makeSyntheticData(mtrue,std)
+    survey.std = std 
+    survey.eps = 1e-5*np.linalg.norm(survey.dobs)
 
     if plotIt:
         import matplotlib.pyplot as plt
         fig, ax = plt.subplots(1,1, figsize = (10, 6))
-        ax.loglog(rx.times, dtrue, 'b.-')
+        ax.loglog(rx.times, survey.dtrue, 'b.-')
         ax.loglog(rx.times, survey.dobs, 'r.-')
         ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
         ax.set_xlabel('Time (s)', fontsize = 14)
@@ -76,6 +74,7 @@ def run(plotIt=True):
     reg = Regularization.Tikhonov(regMesh)
     opt = Optimization.InexactGaussNewton(maxIter = 5)
     invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
     # Create an inversion object
     beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
     betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)

SimPEG/Examples/__init__.py

@@ -1,7 +1,9 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
+import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_FDEM_SusEffects
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
@@ -13,8 +15,9 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
+import MT_1D_analytic_nlayer_Earth
 
-__examples__ = ["EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
+__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Exceptions.py

@@ -1,13 +0,0 @@
-
-
-class SimPEGException(Exception):
-
-    def __init__(self, reason=''):
-        self.reason = reason
-
-    def __str__(self):
-        return '%s: %s' %(self.__class__.__name__,  self.reason)
-
-
-class PairingException(SimPEGException):
-    pass

SimPEG/Maps.py

@@ -4,7 +4,6 @@ from Tests import checkDerivative
 from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
-from scipy.spatial import cKDTree
 
 class IdentityMap(object):
     """
@@ -122,7 +121,7 @@ class IdentityMap(object):
             if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
-        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray! Not a %s'%type(val))
+        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
 
     def __str__(self):
         return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
@@ -459,89 +458,6 @@ class Mesh2Mesh(IdentityMap):
         return self.P
 
 
-class Mesh2MeshTopo(IdentityMap):
-    """
-        Takes a model on one mesh are translates it to another mesh
-        with consideration of topography
-
-    """
-    tree = None
-    nIterpPts = 6
-    P = None #: The CSR projection matrix.
-
-    def __init__(self, meshes, actinds, **kwargs):
-        Utils.setKwargs(self, **kwargs)
-
-        assert type(meshes) is list, "meshes must be a list of two meshes"
-        assert len(meshes) == 2, "meshes must be a list of two meshes"
-        assert type(actinds) is list, "actinds must be a list of two meshes"
-        assert len(actinds) == 2, "actinds must be a list of two meshes"
-        assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
-
-        self.mesh  = meshes[0]
-        self.mesh2 = meshes[1]
-        self.actind = actinds[0]
-        self.actind2 = actinds[1]
-        self._createProjection()
-
-        # Old version using SimPEG interpolation
-        # self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
-
-    def genActiveindfromTopo(mesh, xyztopo):
-        #TODO: This possibly needs to be improved use vtk(?)
-        if mesh.dim==3:
-            nCxy = mesh.nCx*mesh.nCy
-            Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F')
-            Ftopo = NearestNDInterpolator(xyztopo[:,:2], xyztopo[:,2])
-            XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
-            XY.shape
-            topo = Ftopo(XY)
-            actind = []
-            for ixy in range(nCxy):
-                actind.append(topo[ixy] <= Zcc[ixy,:])
-        else:
-            raise NotImplementedError("Only 3D is working")
-
-        return Utils.mkvc(np.vstack(actind))
-
-    #Question .. is it only generated once?
-    def _createProjection(self):
-        """
-            KD Tree interpolation onto the active cells.
-        """
-        if self.tree==None:
-            self.tree = cKDTree(zip(self.mesh.gridCC[self.actind,0], self.mesh.gridCC[self.actind,1], self.mesh.gridCC[self.actind,2]))
-        d, inds = self.tree.query(zip(self.mesh2.gridCC[self.actind2,0],self.mesh2.gridCC[self.actind2,1],self.mesh2.gridCC[self.actind2,2]), k=self.nIterpPts)
-        # Not sure consideration of the volume ...
-        # vol = np.zeros((self.actind2.sum(), self.nIterpPts))
-        # for i in range(self.nIterpPts):
-        #     vol[:,i] = self.mesh.vol[inds[:,i]]
-        w = 1. / d**2
-        w = Utils.sdiag(1./np.sum(w, axis=1)) * (w)
-        I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(self.nIterpPts, axis=1))
-        J = Utils.mkvc(inds)
-        P = sp.coo_matrix( (Utils.mkvc(w),(I, J)), shape=(inds.shape[0], (self.actind).sum()) )
-        # self.P = Utils.sdiag(self.mesh2.vol[self.actind2])*P.tocsc()
-        self.P = P.tocsr()
-
-    @property
-    def shape(self):
-        """Number of parameters in the model."""
-        # return (self.mesh.nC, self.mesh2.nC)
-        return (self.actind2.sum(), self.actind.sum())
-
-    @property
-    def nP(self):
-        """Number of parameters in the model."""
-        # return self.mesh2.nC
-        return self.actind2.sum()
-
-    def _transform(self, m):
-        return self.P*m
-
-    def deriv(self, m):
-        return self.P
-
 class ActiveCells(IdentityMap):
     """
         Active model parameters.

SimPEG/Mesh/DiffOperators.py

@@ -565,58 +565,7 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
-    
-    def unitCellGradx():
-        doc = """Cell centered Gradient in the x dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradx', None) is None:
-                
-                n = self.vnC
-                gx = ddx(n[0]-1)
-                gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
-                
-                self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square)
-        
-            return self._unitCellGradx
-        return locals()
-    unitCellGradx = property(**unitCellGradx())
-    
-    def unitCellGrady():
-        doc = """Cell centered Gradient in they dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGrady', None) is None:
-        
-                n = self.vnC
-                gy = ddx(n[1]-1)
-                gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
-                
-                self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0]))
-            
-            return self._unitCellGrady
-        return locals()
-    unitCellGrady = property(**unitCellGrady())
-        
-    def unitCellGradz():
-        doc = """Cell centered Gradient in they dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradz', None) is None:
-        
-                n = self.vnC
-                gz = ddx(n[2]-1)
-                gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
-                
-                self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0]))
-            
-            return self._unitCellGradz
-        return locals()
-    unitCellGradz = property(**unitCellGradz())        
-        
+
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Optimization.py

@@ -989,19 +989,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
                     cgFlag = 1
                 # End CG Iterations
-                    
+
             # Take a gradient step on the active cells if exist
             if temp != self.xc.size:
-        
+
                 rhs_a = (Active) * -self.g
-                
+
                 dm_i = max( abs( delx ) )
-                dm_a = max( abs(rhs_a) )           
-    
+                dm_a = max( abs(rhs_a) )
+
                 delx = delx + rhs_a * dm_i / dm_a /10.
-                    
+
             # Only keep gradients going in the right direction on the active set
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.
 
-            return delx
+            return delx

SimPEG/Problem.py

@@ -1,6 +1,6 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
-import Maps, Mesh, Exceptions
+import Maps, Mesh
 from Fields import Fields, TimeFields
 
 class BaseProblem(object):
@@ -18,13 +18,9 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    PropMap = None    #: A SimPEG PropertyMap class.
+    mesh    = None    #: A SimPEG.Mesh instance.
 
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Utils.setKwargs(self, **kwargs)
-        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mesh = mesh
-        self.mapping = mapping or Maps.IdentityMap(mesh)
+    PropMap = None    #: A SimPEG PropertyMap class.
 
     @property
     def mapping(self):
@@ -36,9 +32,14 @@ class BaseProblem(object):
             val._assertMatchesPair(self.mapPair)
             self._mapping = val
         else:
-            self._propMapMapping = val
             self._mapping = self.PropMap(val)
 
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.mesh = mesh
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+
     @property
     def survey(self):
         """
@@ -46,22 +47,13 @@ class BaseProblem(object):
         """
         return getattr(self, '_survey', None)
 
-    def pair(self, survey):
+    def pair(self, d):
         """Bind a survey to this problem instance using pointers."""
-        assert isinstance(survey, self.surveyPair), "Survey must be an instance of a %s class."%(self.surveyPair.__name__)
-        if survey.ispaired:
+        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
+        if d.ispaired:
             raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
-        try:
-            self._survey = survey
-            self._validatePairing()
-        except Exceptions.PairingException, e:
-            self._survey = None
-            raise e
-        survey._prob = self
-
-    def _validatePairing(self):
-        """Called when the pair is done, raise a SimPEG.Exceptions.PairingException if unsuccessful"""
-        pass
+        self._survey = d
+        d._prob = self
 
     def unpair(self):
         """Unbind a survey from this problem instance."""
@@ -166,9 +158,6 @@ class BaseProblem(object):
 
 class BaseTimeProblem(BaseProblem):
     """Sets up that basic needs of a time domain problem."""
-    
-    waveformType = "STEPOFF"    
-    current = None
 
     @property
     def timeSteps(self):
@@ -195,11 +184,6 @@ class BaseTimeProblem(BaseProblem):
         self._timeSteps = Utils.meshTensor(value)
         del self.timeMesh
 
-    def currentwaveform(self, wave):
-        self._timeSteps = np.diff(wave[:,0])
-        self.current = wave[:,1]
-        self.waveformType = "GENERAL"
-
     @property
     def nT(self):
         "Number of time steps."
@@ -230,191 +214,4 @@ class BaseTimeProblem(BaseProblem):
             del self._timeMesh
 
 
-class GlobalProblem(BaseProblem):
-    """
-
-        The GlobalProblem allows you to run a whole bunch of SubProblems,
-        potentially in parallel, potentially of different meshes.
-
-        This is handy for working with lots of sources,
-
-    """
-
-    surveyKwargs = {}
-    probKwargs   = {}
-
-    def __init__(self, SubProblem, globalMesh, mapping=None, **kwargs):
-
-        # assert isclass??(SubProblem, BaseProblem), "SubProblem must be a SimPEG.Problem.BaseProblem object."
-        self.surveyPair = SubProblem.surveyPair
-        self.PropMap = SubProblem.PropMap
-        self.mapPair = SubProblem.mapPair
-        self.SubProblem = SubProblem
-
-        Utils.setKwargs(self, **kwargs)
-        assert isinstance(globalMesh, Mesh.BaseMesh), "globalMesh must be a SimPEG.Mesh object."
-        self.globalMesh = globalMesh
-        self.mapping = mapping or Maps.IdentityMap()
-
-    @property
-    def groups(self):
-        """
-            List of lists/integers to say how the sources are grouped.
-
-            e.g.
-
-            survey.srcList = [s0,s1,s2,s3,s4]
-            groups = [ [0,4], [1,3], 2 ]
-        """
-        if getattr(self, '_groups', None) is None:
-            if not self.ispaired: return None
-            self._groups = range(self.survey.nSrc)
-        return self._groups
-    @groups.setter
-    def groups(self, val):
-        assert type(val) is list, 'This should be an list of groups'
-        if self.ispaired:
-            for g in val:
-                assert type(g) in [int, list], 'Must be an integer or a list'
-                if type(g) is int:
-                    assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
-                if type(g) is list:
-                    for sg in g:
-                        assert type(g) is int, 'Must be an integer or a list'
-                        assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
-            assert len(val) == len(self.survey.srcList), 'The groups must be the same length as the srcList in the survey'
-        self._groups = val
-        self._nGroups = None
-
-    @property
-    def meshes(self):
-        if getattr(self, '_meshes', None) is None:
-            if not self.ispaired: return None
-            self._meshes = [self.globalMesh]*self.nGroups
-        return self._meshes
-    @meshes.setter
-    def meshes(self, val):
-        assert type(val) is list
-        if self.ispaired:
-            assert len(val) == self.nGroups
-        self._meshes = val
-
-    @property
-    def nGroups(self):
-        if getattr(self, '_groups', None) is None:
-            return None
-        return len(self.groups)
-
-    def _validatePairing(self):
-        try:
-            self.groups = self.groups # check the assumptions for the grouping
-        except Exception, e:
-            raise Exceptions.PairingException(reason='The grouping does not match the survey')
-        if self.nGroups is not len(self.meshes):
-            raise Exceptions.PairingException(reason='The meshes are not the the same length as the number of groups')
-
-    def getSubProblemandSubSurvey(self, subMap, ind):
-        #This is a core place that we can proceed parallelization
-        assert self.ispaired, 'You must be paired to a survey'
-        assert type(ind) in [int,long] and ind >= 0 and ind < self.nGroups, 'ind must be an index into the group list'
-
-        subMesh = self.meshes[ind]
-        # subMap  = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
-        # subMap  = self.getSubMap(subMesh, ind)
-
-        if self.PropMap is None:
-            prob = self.SubProblem(subMesh, mapping=subMap * self.mapping, **self.probKwargs)
-        else:
-            # This will not work with a fancier propmap...
-            prob = self.SubProblem(subMesh, mapping=subMap * self._propMapMapping, **self.probKwargs)
-
-        survey = self.survey.__class__(srcList=self.survey.srcList[self.groups[ind]], **self.surveyKwargs)
-        prob.pair(survey)
-
-        return prob, survey
-
-    # Not sure we need this here ... 
-
-    def getSubMap(self, subMesh, ind):
-        """The sub"""
-        mesh2mesh = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
-
-        if self.PropMap is None:
-            subMap = mesh2mesh * self.mapping
-        else:
-            subMap = mesh2mesh * self._propMapMapping
-
-        return subMap
-
-if __name__ == '__main__':
-
-
-    from SimPEG import *
-    from SimPEG import EM
-    from scipy.constants import mu_0
-    from pymatsolver import MumpsSolver
-
-    cs = 10.
-    ncx, ncy, ncz = 10, 10, 10
-    npad = 4
-    freq = 1e2
-
-    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
-    mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
-
-    mapping = Maps.ExpMap(mesh)
-
-    x = np.linspace(-10,10,5)
-    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-    rxList = EM.FDEM.Rx(XYZ, 'exi')
-    Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
-    Src1 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
-
-
-    prb0 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
-    survey = EM.FDEM.Survey([Src0])
-    prb0.pair(survey)
-    prb1 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
-    survey = EM.FDEM.Survey([Src1])
-    prb1.pair(survey)
-
-
-
-    sig = 1e-1
-    sigma = np.ones(mesh.nC)*sig
-    sigma[mesh.gridCC[:,2] > 0] = 1e-8
-    m = np.log(sigma)
-
-    GP = GlobalProblem(EM.FDEM.Problem_b, mesh, mapping=mapping, meshes=[mesh,mesh])
-    survey = EM.FDEM.Survey([Src0, Src1])
-    GP.pair(survey)
-
-    gp1 = GP.getSubProblem(0)
-    gp1.Solver = MumpsSolver
-
-    pu = prb0.fields(m)
-    gpu = gp1.fields(m)
-
-    bfz = mesh.r(pu[Src0, 'b'],'F','Fz','M')
-    bfz = mesh.r(gpu[Src0, 'b'],'F','Fz','M')
-    x = np.linspace(-55,55,12)
-    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-    P = mesh.getInterpolationMat(XYZ, 'Fz')
-
-    # an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, Src0.freq, sig)
-
-    # diff = np.log10(np.abs(P*np.imag(pu[Src0, 'b']) - mu_0*np.imag(an)))
-    # diff = np.log10(np.abs(P*np.imag(gpu[Src0, 'b']) - mu_0*np.imag(an)))
-
-    import matplotlib.pyplot as plt
-    plt.plot(x,np.log10(np.abs(P*np.imag(pu[Src0, 'b']))), 'r-s')
-    plt.plot(x,np.log10(np.abs(P*np.imag(gpu[Src0, 'b']))), 'b')
-    # plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
-    # plt.plot(x,diff,'g')
-    plt.show()
-
-
-
 

SimPEG/Regularization.py

@@ -282,239 +282,3 @@ class Tikhonov(BaseRegularization):
             out = mD.T * ( self.W.T * r )
         return out
 
-class Simple(BaseRegularization):
-    """
-        Only for tensor mesh
-    """     
-    
-    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
-    alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
-    alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
-    alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
-
-    def __init__(self, mesh, mapping=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
-
-    
-        
-    @property
-    def Ws(self):
-        """Regularization matrix Ws"""
-        if getattr(self,'_Ws', None) is None:
-                self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-        return self._Ws
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
-        return self._Wz
-
-    @property
-    def Wxx(self):
-        """Regularization matrix Wxx"""
-        if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-        return self._Wxx
-
-    @property
-    def Wyy(self):
-        """Regularization matrix Wyy"""
-        if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-        return self._Wyy
-
-    @property
-    def Wzz(self):
-        """Regularization matrix Wzz"""
-        if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-        return self._Wzz
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def eval(self, m):
-        if self.smoothModel == True:
-            r1 = self.Wsmooth * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (m - self.mref) )
-            return 0.5*(r1.dot(r1)+r2.dot(r2))
-        elif self.smoothModel == False:
-            r = self.W * ( self.mapping * (m - self.mref) )
-            return 0.5*r.dot(r)
-
-
-    @Utils.timeIt
-    def evalDeriv(self, m):
-        """
-
-        The regularization is:
-
-        .. math::
-
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
-
-        So the derivative is straight forward:
-
-        .. math::
-
-            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
-        """
-        if self.smoothModel == True:
-            mD1 = self.mapping.deriv(m)
-            mD2 = self.mapping.deriv(m - self.mref)
-            r1 = self.Wsmooth * ( self.mapping * (m))
-            r2 = self.Ws * ( self.mapping * (m - self.mref) )
-            out1 = mD1.T * ( self.Wsmooth.T * r1 )
-            out2 = mD2.T * ( self.Ws.T * r2 )
-            out = out1+out2
-        elif self.smoothModel == False:
-            mD = self.mapping.deriv(m - self.mref)
-            r = self.W * ( self.mapping * (m - self.mref) )
-            out = mD.T * ( self.W.T * r )
-        return out
-
-class SparseRegularization(Simple):
-    
-    eps = 1e-1
-
-    m = None
-    gamma = 1.
-    p = 0.
-    qx = 2.
-    qy = 2.
-    qz = 2.
-    
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, **kwargs)
-        
-    
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-    
-    @property
-    def Ws(self):
-        """Regularization matrix Ws"""
-        if getattr(self, 'm', None) is None:
-            self.Rs = Utils.speye(self.mesh.nC)            
-            
-        else:
-            f_m = self.m
-            self.rs = self.R(f_m , self.p, self.eps)
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-            
-        self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
-            
-        return self._Ws
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        
-        if getattr(self, 'm', None) is None:
-            self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGradx * self.m
-            self.rx = self.R( f_m , self.qx, self.eps)
-            self.Rx = Utils.sdiag( self.rx )
-            
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        
-        if getattr(self, 'm', None) is None:
-            self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGrady * self.m
-            self.ry = self.R( f_m , self.qy, self.eps)
-            self.Ry = Utils.sdiag( self.ry )
-            
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'm', None) is None:
-            self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGradz * self.m
-            self.rz = self.R( f_m , self.qz, self.eps)
-            self.Rz = Utils.sdiag( self.rz )
-        
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
-        return self._Wz    
-
-    
-    def R(self, f_m , p, dec):
-
-        eta = (self.eps**(1-p/2.))**0.5
-        r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
-
-        return r

SimPEG/Survey.py

@@ -1,5 +1,6 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
 
+
 class BaseRx(object):
     """SimPEG Receiver Object"""
 
@@ -204,6 +205,7 @@ class BaseSurvey(object):
     __metaclass__ = Utils.SimPEGMetaClass
 
     std = None       #: Estimated Standard Deviations
+    eps = None       #: Estimated Noise Floor
     dobs = None      #: Observed data
     dtrue = None     #: True data, if data is synthetic
     mtrue = None     #: True model, if data is synthetic
@@ -222,8 +224,6 @@ class BaseSurvey(object):
 
     @srcList.setter
     def srcList(self, value):
-        if isinstance(value, self.srcPair):
-            value = [value]
         assert type(value) is list, 'srcList must be a list'
         assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
         assert len(set(value)) == len(value), 'The srcList must be unique'
@@ -375,11 +375,3 @@ class BaseSurvey(object):
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
-
-class LinearSurvey(BaseSurvey):
-    def projectFields(self, u):
-        return u
-    
-    @property
-    def nD(self):
-        return self.prob.G.shape[0]

SimPEG/Utils/ModelBuilder.py

@@ -118,44 +118,6 @@ def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.
     D = np.sqrt(np.sum(G**2,axis=1))
     return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5
 
-def getIndicesSphere(center,radius,ccMesh):
-    """
-        Creates a vector containing the sphere indices in the cell centers mesh.
-        Returns a tuple
-
-        The sphere is defined by the points
-
-        p0, describe the position of the center of the cell
-
-        r, describe the radius of the sphere.
-
-        ccMesh represents the cell-centered mesh
-
-        The points p0 must live in the the same dimensional space as the mesh.
-
-    """
-
-    # Validation: mesh and point (p0) live in the same dimensional space
-    dimMesh = np.size(ccMesh[0,:])
-    assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
-
-    if dimMesh == 1:
-       # Define the reference points
-        
-        ind  = np.abs(center[0] - ccMesh[:,0]) < radius
-
-    elif dimMesh == 2:
-       # Define the reference points
-
-        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius
-
-    elif dimMesh == 3:
-        # Define the points
-        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius
-
-    # Return a tuple
-    return ind
-
 def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     """
     Define a two layered model.  Depth of the first layer must be specified.

SimPEG/Utils/matutils.py

@@ -2,7 +2,6 @@ import numpy as np
 import scipy.sparse as sp
 from codeutils import isScalar
 
-
 def mkvc(x, numDims=1):
     """Creates a vector with the number of dimension specified
 
@@ -26,6 +25,9 @@ def mkvc(x, numDims=1):
     if hasattr(x, 'tovec'):
         x = x.tovec()
 
+    if isinstance(x, Zero):
+        return x
+    
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -37,6 +39,9 @@ def mkvc(x, numDims=1):
 
 def sdiag(h):
     """Sparse diagonal matrix"""
+    if isinstance(h, Zero):
+        return Zero()
+
     return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
 
 def sdInv(M):
@@ -417,6 +422,12 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
+    @property 
+    def transpose(self): return Zero()
+    
+    @property
+    def T(self): return Zero()
+
 class Identity(object):
     _positive = True
     def __init__(self, positive=True):

docs/em/api_FDEM.rst

@@ -19,14 +19,14 @@ Electromagnetic phenomena are governed by Maxwell's equations. They describe the
 
 Fourier Transform Convention
 ----------------------------
-In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the \\(e^{i \\omega t} \\) convention, so we define our Fourier Transform pair as
+In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the :math:`e^{i \omega t}` convention, so we define our Fourier Transform pair as
 
 .. math ::
-	F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\
+    F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\
 
-	f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega
+    f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega
 
-where \\(\\omega\\) is angular frequency, \\(t\\) is time, \\(F(\\omega)\\) is the function defined in the frequency domain and \\(f(t)\\) is the function defined in the time domain.
+where :math:`\omega` is angular frequency, :math:`t` is time, :math:`F(\omega)` is the function defined in the frequency domain and :math:`f(t)` is the function defined in the time domain.
 
 
 Maxwell's Equations
@@ -34,44 +34,46 @@ Maxwell's Equations
 In the frequency domain, Maxwell's equations are given by
 
 .. math ::
-	\curl \vec{E} = - i \omega \vec{B} \\
+    \curl \vec{E} + i \omega \vec{B} = \vec{S_m}\\
 
-	\curl \vec{H} = \vec{J} + i \omega \vec{D} + \vec{S} \\
+    \curl \vec{H} - \vec{J} - i \omega \vec{D} = \vec{S_e} \\
 
-	\div \vec{B} = 0 \\
+    \div \vec{B} = 0 \\
 
-	\div \vec{D} = \rho_f
+    \div \vec{D} = \rho_f
 
 where:
 
-	- \\(\\vec{E}\\) : electric field (\\(V/m\\))
-	- \\(\\vec{H}\\) : magnetic field (\\(A/m\\))
-	- \\(\\vec{B}\\) : magnetic flux density (\\(Wb/m^2\\))
-	- \\(\\vec{D}\\) : electric displacement / electric flux density (\\(C/m^2\\))
-	- \\(\\vec{J}\\) : electric current density (\\(A/m^2\\))
-	- \\(\\rho_f\\) : free charge density
+    - :math:`\vec{E}`   : electric field (:math:`V/m` )
+    - :math:`\vec{H}`   : magnetic field (:math:`A/m` )
+    - :math:`\vec{B}`   : magnetic flux density (:math:`Wb/m^2` )
+    - :math:`\vec{D}`   : electric displacement / electric flux density (:math:`C/m^2` )
+    - :math:`\vec{J}`   : electric current density (:math:`A/m^2` )
+    - :math:`\vec{S_m}` : magnetic source term (:math:`V/m^2` )
+    - :math:`\vec{S_e}` : electric source term (:math:`A/m^2` )
+    - :math:`\rho_f`   : free charge density (:math:`\Omega m` )
 
-The source term is \\(\\vec{S}\\)
 
 
 Constitutive Relations
 ----------------------
+
 The fields and fluxes are related through the constitutive relations. At each frequency, they are given by
 
 .. math ::
-	\vec{J} = \sigma \vec{E} \\
+    \vec{J} = \sigma \vec{E} \\
 
-	\vec{B} = \mu \vec{H} \\
+    \vec{B} = \mu \vec{H} \\
 
-	\vec{D} = \varepsilon \vec{E}
+    \vec{D} = \varepsilon \vec{E}
 
 where:
 
-- \\(\\sigma\\) : electrical conductivity \\(S/m\\)
-- \\(\\mu\\) : magnetic permeability \\(H/m\\)
-- \\(\\varepsilon\\) : dielectric permittivity \\(F/m\\)
+- :math:`\sigma` : electrical conductivity (:math:`S/m`)
+- :math:`\mu` : magnetic permeability (:math:`H/m`)
+- :math:`\varepsilon` : dielectric permittivity (:math:`F/m`)
 
-\\(\\sigma\\), \\(\\mu\\), \\(\\varepsilon\\) are physical properties which depend on the material. \\(\\sigma\\) describes how easily electric current passes through a material, \\(\\mu\\) describes how easily a material is magnetized, and \\(\\varepsilon\\) describes how easily a material is electrically polarized. In most geophysical applications of EM, \\(\\sigma\\) is the the primary physical property of interest, and \\(\\mu\\), \\(\\varepsilon\\) are assumed to have their free-space values \\(\\mu_0 = 4\\pi \\times 10^{-7} H/m \\), \\(\\varepsilon_0 = 8.85 \\times 10^{-12} F/m\\)
+:math:`\sigma`, :math:`\mu`, :math:`\varepsilon` are physical properties which depend on the material. :math:`\sigma` describes how easily electric current passes through a material, :math:`\mu` describes how easily a material is magnetized, and :math:`\varepsilon` describes how easily a material is electrically polarized. In most geophysical applications of EM, :math:`\sigma` is the the primary physical property of interest, and :math:`\mu`, :math:`\varepsilon` are assumed to have their free-space values :math:`\mu_0 = 4\pi \times 10^{-7} H/m` , :math:`\varepsilon_0 = 8.85 \times 10^{-12} F/m`
 
 
 Quasi-static Approximation
@@ -80,8 +82,8 @@ Quasi-static Approximation
 For the frequency range typical of most geophysical surveys, the contribution of the electric displacement is negligible compared to the electric current density. In this case, we use the Quasi-static approximation and assume that this term can be neglected, giving
 
 .. math ::
-	\nabla \times \vec{E} = -i \omega \vec{B} \\
-	\nabla \times \vec{H} = \vec{J} + \vec{S}
+    \nabla \times \vec{E} + i \omega \vec{B} = \vec{S_m} \\
+    \nabla \times \vec{H} - \vec{J} = \vec{S_e}
 
 
 Implementation in SimPEG.EM
@@ -90,14 +92,14 @@ Implementation in SimPEG.EM
 We consider two formulations in SimPEG.EM, both first-order and both in terms of one field and one flux. We allow for the definition of magnetic and electric sources (see for example: Ward and Hohmann, starting on page 144). The E-B formulation is in terms of the electric field and the magnetic flux:
 
 .. math ::
-	\nabla \times \vec{E} + i \omega \vec{B} = \vec{S}_m \\
-	\nabla \times \mu^{-1} \vec{B} - \sigma \vec{E} = \vec{S}_e
+    \nabla \times \vec{E} + i \omega \vec{B} = \vec{S}_m \\
+    \nabla \times \mu^{-1} \vec{B} - \sigma \vec{E} = \vec{S}_e
 
 The H-J formulation is in terms of the current density and the magnetic field:
 
 .. math ::
-	\nabla \times \sigma^{-1} \vec{J} + i \omega \mu \vec{H} = \vec{S}_m \\
-	\nabla \times \vec{H} - \vec{J} = \vec{S}_e
+    \nabla \times \sigma^{-1} \vec{J} + i \omega \mu \vec{H} = \vec{S}_m \\
+    \nabla \times \vec{H} - \vec{J} = \vec{S}_e
 
 
 Discretizing
@@ -106,34 +108,34 @@ For both formulations, we use a finite volume discretization
 and discretize fields on cell edges, fluxes on cell faces and
 physical properties in cell centers. This is particularly
 important when using symmetry to reduce the dimensionality of a problem
-(for instance on a 2D CylMesh, there are \\(r\\), \\(z\\) faces and \\(\\theta\\) edges)
+(for instance on a 2D CylMesh, there are :math:`r`, :math:`z` faces and :math:`\theta` edges)
 
 .. figure:: ../images/finitevolrealestate.png
-	:align: center
-	:scale: 60 %
+    :align: center
+    :scale: 60 %
 
 For the two formulations, the discretization of the physical properties, fields and fluxes are summarized below.
 
 .. figure:: ../images/ebjhdiscretizations.png
-	:align: center
-	:scale: 60 %
+    :align: center
+    :scale: 60 %
 
-Note that resistivity is the inverse of conductivity, \\(\\rho = \\sigma^{-1}\\).
+Note that resistivity is the inverse of conductivity, :math:`\rho = \sigma^{-1}`.
 
 
-E-B Formulation:
-****************
+E-B Formulation
+---------------
 
 .. math ::
-	\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\
-	\mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}
+    \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\
+    \mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}
 
-H-J Formulation:
-****************
+H-J Formulation
+---------------
 
 .. math ::
-	\mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\
-	\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
+    \mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\
+    \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
 
 .. Forward Problem
@@ -144,6 +146,10 @@ H-J Formulation:
 
 API
 ===
+
+FDEM Problem
+------------
+
 .. automodule:: SimPEG.EM.FDEM.FDEM
     :show-inheritance:
     :members:
@@ -157,3 +163,17 @@ FDEM Survey
     :show-inheritance:
     :members:
     :undoc-members:
+
+.. automodule:: SimPEG.EM.FDEM.SrcFDEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+FDEM Fields
+-----------
+
+.. automodule:: SimPEG.EM.FDEM.FieldsFDEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+

docs/em/api_TDEM.rst

@@ -48,6 +48,305 @@
     \newcommand{\I}{\vec{I}}
 
 
+Time Domain Electromagnetics
+****************************
+
+.. _api_TDEM_derivation:
+
+Time-Domain EM Derivation
+=========================
+
+The following shows the derivation for the TDEM problem. We use the b-formulation below.
+(More to come soon..!)
+
+
+Sensitivity Calculation
+-----------------------
+
+.. math::
+
+    \begin{align}
+        \dcurl \e^{(t+1)} + \frac{\b^{(t+1)} - \b^{(t)}}{\delta t} = 0 \\
+        \dcurl^\top \MfMui \b^{(t+1)} - \MeSig \e^{(t+1)} = \Me \j_s^{(t+1)}
+    \end{align}
+
+Using Gauss-Newton to solve the inverse problem requires the ability to calculate the product of the
+Jacobian and a vector, as well as the transpose of the Jacobian times a vector.
+The above system can be rewritten as:
+
+.. math::
+
+    \begin{align}
+        \mathbf{A} \u^{(t+1)} + \mathbf{B} \u^{(t)}= \s^{(t+1)}
+    \end{align}
+
+where
+
+.. math::
+
+    \begin{align}
+        \mathbf{A} =
+        \left[
+            \begin{array}{cc}
+                \frac{1}{\delta t} \MfMui & \MfMui\dcurl \\
+                \dcurl^\top \MfMui & -\MeSig
+            \end{array}
+        \right] \\
+        \mathbf{B} =
+        \left[
+            \begin{array}{cc}
+                -\frac{1}{\delta t} \MfMui & 0 \\
+                0 & 0
+            \end{array}
+        \right] \\
+        \u^{(k)} = \left[
+        \begin{array}{c}
+            \b^{(k)}\\
+            \e^{(k)}
+        \end{array}
+        \right] \\
+        \s^{(k)} = \left[
+        \begin{array}{c}
+            0\\
+            \Me \j^{(k)}_s
+        \end{array}
+        \right]
+    \end{align}
+
+.. note::
+
+    Here we have multiplied through by \\(\\MfMui\\) to make A and B symmetric!
+
+The entire time dependent system can be written in a single matrix expression
+
+.. math::
+
+    \begin{align}
+        \hat{\mathbf{A}} \hat{u} = \hat{s}
+    \end{align}
+
+where
+
+.. math::
+
+    \begin{align}
+        \mathbf{\hat{A}} = \left[
+        \begin{array}{cccc}
+            A & 0 & & \\
+            B & A & & \\
+              & \ddots & \ddots & \\
+              & & B & A
+        \end{array}
+        \right] \\
+        \hat{u} = \left[
+            \begin{array}{c}
+                \u^{(1)} \\
+                \u^{(2)} \\
+                \vdots \\
+                \u^{(N)}
+            \end{array} \right]\\
+        \hat{s} = \left[
+            \begin{array}{c}
+                \s^{(1)} - \mathbf{B} \u^{(0)} \\
+                \s^{(2)} \\
+                \vdots \\
+                \s^{(N)}
+            \end{array}
+        \right]
+    \end{align}
+
+For the fields \\(\\u\\), the measured data is given by
+
+.. math::
+
+    \begin{align}
+        \vec{d} = \mathbf{Q} \u
+    \end{align}
+
+The sensitivity matrix **J** is then defined as
+
+.. math::
+
+    \begin{align}
+        \mathbf{J} = \mathbf{Q} \frac{\partial \u}{\partial \sigma}
+    \end{align}
+
+
+Defining the function \\(\\c(m,\\u)\\) to be
+
+.. math::
+
+    \begin{align}
+        \vec{c}(m,\u) = \hat{\mathbf{A}} \vec{u} - \vec{q} = \vec{0}
+    \end{align}
+
+then
+
+.. math::
+
+    \begin{align}
+        \frac{\partial \vec{c}}{\partial m} \partial m
+        + \frac{\partial \vec{c}}{\partial \u} \partial \vec{u} = 0
+    \end{align}
+
+or
+
+.. math::
+
+    \begin{align}
+        \frac{\partial \vec{u}}{\partial m} = -\left(\frac{\partial \vec{c}}{\partial \u} \right)^{-1} \frac{\partial \vec{c}}{\partial m}
+    \end{align}
+
+
+Differentiating, we find that
+
+.. math::
+
+    \begin{align}
+        \frac{\partial \vec{c}}{\partial \hat{u}} = \hat{\mathbf{A}}
+    \end{align}
+
+and
+
+.. math::
+
+    \begin{align}
+        \frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma =
+        \left[
+            \begin{array}{c}
+                g_\sigma^{(1)}\\
+                g_\sigma^{(2)}\\
+                \vdots \\
+                g_\sigma^{(N)}
+            \end{array}
+        \right]
+    \end{align}
+
+with
+
+.. math::
+
+    \begin{align}
+        g_\sigma^{(n)} =
+        \left[
+            \begin{array}{c}
+                \mathbf{0} \\
+                - \diag{\e^{(n)}} \Ace \diag{\vec{V}}
+            \end{array}
+        \right]
+    \end{align}
+
+
+Implementing **J** times a vector
+---------------------------------
+
+Multiplying **J** onto a vector can be broken into three steps
+
+
+* Compute \\(\\vec{p} = \\mathbf{G}m\\)
+* Solve \\(\\hat{\\mathbf{A}} \\vec{y} = \\vec{p}\\)
+* Compute \\(\\vec{w} = -\\mathbf{Q} \\vec{y}\\)
+
+.. math::
+
+    \begin{align}
+        \vec{p}^{(n)} = \left[
+            \begin{array}{c}
+                \vec{p}_b^{(n)} \\
+                \vec{p}_e^{(n)}
+            \end{array}
+        \right] \\
+        \vec{p}_b^{(n)} = 0 \\
+        \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m
+    \end{align}
+
+
+For all time steps:
+
+.. math::
+
+    \begin{align}
+        \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t+1)} + \MfMui\dcurl \vec{y}_{e}^{(t+1)}
+        - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)}
+        = \vec{p}_b^{(t+1)} \\
+        \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig \vec{y}_e^{(t+1)} = \vec{p}_e^{(t+1)}
+    \end{align}
+
+and
+
+.. math::
+
+    \begin{align}
+        \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t+1)} =
+        \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)}
+        + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} + \vec{p}_b^{(t+1)} \\
+        \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
+    \end{align}
+
+.. note::
+
+    For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero.
+
+
+
+
+Implementing **J** transpose times a vector
+-------------------------------------------
+
+Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps
+
+
+* Compute \\(\\vec{p} = \\mathbf{Q}^\\top \\vec{v}\\)
+* Solve \\(\\hat{\\mathbf{A}}^\\top \\vec{y} = \\vec{p}\\)
+* Compute \\(\\vec{w} = -\\mathbf{G}^\\top y\\)
+
+
+.. math::
+
+    \mathbf{\hat{A}}^\top = \left[
+        \begin{array}{cccc}
+            A & B & & \\
+              & \ddots & \ddots & \\
+              & & A & B \\
+              & & 0 & A
+        \end{array}
+    \right]
+
+For the all time-steps (going backwards in time):
+
+
+.. math::
+
+    A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)}
+
+
+.. math::
+
+    \begin{align}
+        \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)}
+        - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)}
+        = \vec{p}_b^{(t)} \\
+        \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)}
+    \end{align}
+
+and
+
+.. math::
+
+    \begin{align}
+        \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} =
+        \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)}
+        + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\
+        \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)}
+    \end{align}
+
+
+.. note::
+
+    For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero.
+
+
+
 TDEM - B formulation
 ====================
 

docs/em/api_TDEM_derivation.rst

@@ -1,341 +0,0 @@
-.. _api_TDEM_derivation:
-
-
-.. math::
-
-    \renewcommand{\div}{\nabla\cdot\,}
-    \newcommand{\grad}{\vec \nabla}
-    \newcommand{\curl}{{\vec \nabla}\times\,}
-    \newcommand {\J}{{\vec J}}
-    \renewcommand{\H}{{\vec H}}
-    \newcommand {\E}{{\vec E}}
-    \newcommand{\dcurl}{{\mathbf C}}
-    \newcommand{\dgrad}{{\mathbf G}}
-    \newcommand{\Acf}{{\mathbf A_c^f}}
-    \newcommand{\Ace}{{\mathbf A_c^e}}
-    \renewcommand{\S}{{\mathbf \Sigma}}
-    \newcommand{\St}{{\mathbf \Sigma_\tau}}
-    \newcommand{\T}{{\mathbf T}}
-    \newcommand{\Tt}{{\mathbf T_\tau}}
-    \newcommand{\diag}[1]{\,{\sf diag}\left( #1 \right)}
-    \newcommand{\M}{{\mathbf M}}
-    \newcommand{\MfMui}{{\M^f_{\mu^{-1}}}}
-    \newcommand{\MeSig}{{\M^e_\sigma}}
-    \newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}}
-    \newcommand{\MeSigO}{{\M^e_{\sigma_0}}}
-    \newcommand{\Me}{{\M^e}}
-    \newcommand{\Mes}[1]{{\M^e_{#1}}}
-    \newcommand{\Mee}{{\M^e_e}}
-    \newcommand{\Mej}{{\M^e_j}}
-    \newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)}
-    \newcommand{\bE}{\mathbf{E}}
-    \newcommand{\bH}{\mathbf{H}}
-    \newcommand{\B}{\vec{B}}
-    \newcommand{\D}{\vec{D}}
-    \renewcommand{\H}{\vec{H}}
-    \newcommand{\s}{\vec{s}}
-    \newcommand{\bfJ}{\bf{J}}
-    \newcommand{\vecm}{\vec m}
-    \renewcommand{\Re}{\mathsf{Re}}
-    \renewcommand{\Im}{\mathsf{Im}}
-    \renewcommand {\j}  { {\vec j} }
-    \newcommand {\h}  { {\vec h} }
-    \renewcommand {\b}  { {\vec b} }
-    \newcommand {\e}  { {\vec e} }
-    \newcommand {\c}  { {\vec c} }
-    \renewcommand {\d}  { {\vec d} }
-    \renewcommand {\u}  { {\vec u} }
-    \newcommand{\I}{\vec{I}}
-
-
-Time-Domain EM Derivation
-*************************
-
-The following shows the derivation for the TDEM problem. We use the b-formulation below.
-(More to come soon..!)
-
-
-Sensitivity Calculation
-=======================
-
-.. math::
-
-    \begin{align}
-        \dcurl \e^{(t+1)} + \frac{\b^{(t+1)} - \b^{(t)}}{\delta t} = 0 \\
-        \dcurl^\top \MfMui \b^{(t+1)} - \MeSig \e^{(t+1)} = \Me \j_s^{(t+1)}
-    \end{align}
-
-Using Gauss-Newton to solve the inverse problem requires the ability to calculate the product of the
-Jacobian and a vector, as well as the transpose of the Jacobian times a vector.
-The above system can be rewritten as:
-
-.. math::
-
-    \begin{align}
-        \mathbf{A} \u^{(t+1)} + \mathbf{B} \u^{(t)}= \s^{(t+1)}
-    \end{align}
-
-where
-
-.. math::
-
-    \begin{align}
-        \mathbf{A} =
-        \left[
-            \begin{array}{cc}
-                \frac{1}{\delta t} \MfMui & \MfMui\dcurl \\
-                \dcurl^\top \MfMui & -\MeSig
-            \end{array}
-        \right] \\
-        \mathbf{B} =
-        \left[
-            \begin{array}{cc}
-                -\frac{1}{\delta t} \MfMui & 0 \\
-                0 & 0
-            \end{array}
-        \right] \\
-        \u^{(k)} = \left[
-        \begin{array}{c}
-            \b^{(k)}\\
-            \e^{(k)}
-        \end{array}
-        \right] \\
-        \s^{(k)} = \left[
-        \begin{array}{c}
-            0\\
-            \Me \j^{(k)}_s
-        \end{array}
-        \right]
-    \end{align}
-
-.. note::
-
-    Here we have multiplied through by \\(\\MfMui\\) to make A and B symmetric!
-
-The entire time dependent system can be written in a single matrix expression
-
-.. math::
-
-    \begin{align}
-        \hat{\mathbf{A}} \hat{u} = \hat{s}
-    \end{align}
-
-where
-
-.. math::
-
-    \begin{align}
-        \mathbf{\hat{A}} = \left[
-        \begin{array}{cccc}
-            A & 0 & & \\
-            B & A & & \\
-              & \ddots & \ddots & \\
-              & & B & A
-        \end{array}
-        \right] \\
-        \hat{u} = \left[
-            \begin{array}{c}
-                \u^{(1)} \\
-                \u^{(2)} \\
-                \vdots \\
-                \u^{(N)}
-            \end{array} \right]\\
-        \hat{s} = \left[
-            \begin{array}{c}
-                \s^{(1)} - \mathbf{B} \u^{(0)} \\
-                \s^{(2)} \\
-                \vdots \\
-                \s^{(N)}
-            \end{array}
-        \right]
-    \end{align}
-
-For the fields \\(\\u\\), the measured data is given by
-
-.. math::
-
-    \begin{align}
-        \vec{d} = \mathbf{Q} \u
-    \end{align}
-
-The sensitivity matrix **J** is then defined as
-
-.. math::
-
-    \begin{align}
-        \mathbf{J} = \mathbf{Q} \frac{\partial \u}{\partial \sigma}
-    \end{align}
-
-
-Defining the function \\(\\c(m,\\u)\\) to be
-
-.. math::
-
-    \begin{align}
-        \vec{c}(m,\u) = \hat{\mathbf{A}} \vec{u} - \vec{q} = \vec{0}
-    \end{align}
-
-then
-
-.. math::
-
-    \begin{align}
-        \frac{\partial \vec{c}}{\partial m} \partial m
-        + \frac{\partial \vec{c}}{\partial \u} \partial \vec{u} = 0
-    \end{align}
-
-or
-
-.. math::
-
-    \begin{align}
-        \frac{\partial \vec{u}}{\partial m} = -\left(\frac{\partial \vec{c}}{\partial \u} \right)^{-1} \frac{\partial \vec{c}}{\partial m}
-    \end{align}
-
-
-Differentiating, we find that
-
-.. math::
-
-    \begin{align}
-        \frac{\partial \vec{c}}{\partial \hat{u}} = \hat{\mathbf{A}}
-    \end{align}
-
-and
-
-.. math::
-
-    \begin{align}
-        \frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma =
-        \left[
-            \begin{array}{c}
-                g_\sigma^{(1)}\\
-                g_\sigma^{(2)}\\
-                \vdots \\
-                g_\sigma^{(N)}
-            \end{array}
-        \right]
-    \end{align}
-
-with
-
-.. math::
-
-    \begin{align}
-        g_\sigma^{(n)} =
-        \left[
-            \begin{array}{c}
-                \mathbf{0} \\
-                - \diag{\e^{(n)}} \Ace \diag{\vec{V}}
-            \end{array}
-        \right]
-    \end{align}
-
-
-Implementing **J** times a vector
-=================================
-
-Multiplying **J** onto a vector can be broken into three steps
-
-
-* Compute \\(\\vec{p} = \\mathbf{G}m\\)
-* Solve \\(\\hat{\\mathbf{A}} \\vec{y} = \\vec{p}\\)
-* Compute \\(\\vec{w} = -\\mathbf{Q} \\vec{y}\\)
-
-.. math::
-
-    \begin{align}
-        \vec{p}^{(n)} = \left[
-            \begin{array}{c}
-                \vec{p}_b^{(n)} \\
-                \vec{p}_e^{(n)}
-            \end{array}
-        \right] \\
-        \vec{p}_b^{(n)} = 0 \\
-        \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m
-    \end{align}
-
-
-For all time steps:
-
-.. math::
-
-    \begin{align}
-        \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t+1)} + \MfMui\dcurl \vec{y}_{e}^{(t+1)}
-        - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)}
-        = \vec{p}_b^{(t+1)} \\
-        \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig \vec{y}_e^{(t+1)} = \vec{p}_e^{(t+1)}
-    \end{align}
-
-and
-
-.. math::
-
-    \begin{align}
-        \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t+1)} =
-        \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)}
-        + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} + \vec{p}_b^{(t+1)} \\
-        \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
-    \end{align}
-
-.. note::
-
-    For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero.
-
-
-
-
-Implementing **J** transpose times a vector
-===========================================
-
-Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps
-
-
-* Compute \\(\\vec{p} = \\mathbf{Q}^\\top \\vec{v}\\)
-* Solve \\(\\hat{\\mathbf{A}}^\\top \\vec{y} = \\vec{p}\\)
-* Compute \\(\\vec{w} = -\\mathbf{G}^\\top y\\)
-
-
-.. math::
-
-    \mathbf{\hat{A}}^\top = \left[
-        \begin{array}{cccc}
-            A & B & & \\
-              & \ddots & \ddots & \\
-              & & A & B \\
-              & & 0 & A
-        \end{array}
-    \right]
-
-For the all time-steps (going backwards in time):
-
-
-.. math::
-
-    A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)}
-
-
-.. math::
-
-    \begin{align}
-        \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)}
-        - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)}
-        = \vec{p}_b^{(t)} \\
-        \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)}
-    \end{align}
-
-and
-
-.. math::
-
-    \begin{align}
-        \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} =
-        \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)}
-        + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\
-        \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)}
-    \end{align}
-
-
-.. note::
-
-    For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero.

docs/em/api_Utils.rst

@@ -4,6 +4,16 @@ simpegEM Utilities
 SimPEG for EM provides a few EM specific utility codes,
 sources, and analytic functions.
 
+Utilities for Electromagnetics
+==============================
+
+.. automodule:: SimPEG.EM.Utils
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
 Analytic Functions - Time
 =========================
 
@@ -22,12 +32,3 @@ Analytic Functions - Frequency
     :members:
     :undoc-members:
     :inherited-members:
-
-
-Sources
-=======
-
-.. autoclass:: SimPEG.EM.FDEM.SrcFDEM.MagDipole
-    :show-inheritance:
-    :members:
-    :undoc-members:

docs/em/index.rst

@@ -3,42 +3,24 @@ Electromagnetics
 ================
 
 `SimPEG.EM` uses SimPEG as the framework for the forward and inverse
-electromagnetics geophysical problems.
+electromagnetics geophysical problems. 
 
-Time Domian Electromagnetics
-----------------------------
-
-.. toctree::
-    :maxdepth: 2
-
-    api_TDEM_derivation
+To solve for predicted data, we follow the framework shown below. The model is
+what we invert for. This is mapped to a physical property on the simulation
+mesh. A source which is used to excite the system is specified. Having a model
+and a source, we can solve Maxwell's equations for fields. We sample these
+fields with recievers to give us predicted data. 
 
 
-Code for Time Domian Electromagnetics
--------------------------------------
+.. image:: ../images/simpegEM_noMath.png
+    :scale: 50%
 
-.. toctree::
-    :maxdepth: 2
-
-    api_TDEM
-
-Frequency Domian Electromagnetics
----------------------------------
 
 .. toctree::
    :maxdepth: 2
 
    api_FDEM
-
-
-Utility Codes
--------------
-
-.. toctree::
-   :maxdepth: 2
-
+   api_TDEM
    api_Utils
 
 
-
-

docs/examples/EM_FDEM_1D_Inversion.rst

@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_1D_Inversion:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: FDEM: 1D: Inversion
+=======================
+
+Here we will create and run a FDEM 1D inversion.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_1D_Inversion.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+    :language: python
+    :linenos:

docs/examples/EM_FDEM_SusEffects.rst

@@ -0,0 +1,41 @@
+.. _examples_EM_FDEM_SusEffects:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: FDEM: Effects of susceptibility
+===================================
+
+When airborne freqeuncy domain EM (AFEM) survey is flown over
+the earth including significantly susceptible bodies (magnetite-rich rocks),
+negative data is often observed in the real part of the lowest frequency
+(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
+occurs due to a susceptible body when the magnetic field is applied.
+
+To clarify what is happening in the earth when we are exciting the earth with
+a loop source in the frequency domain we run three forward modelling:
+
+    - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
+    - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
+    - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
+
+We plot vector magnetic fields in the earth. For secondary fields we provide
+F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
+figure show both real and parts.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_SusEffects.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py
+    :language: python
+    :linenos:

tests/examples/test_examples.py

@@ -1,8 +1,28 @@
 import unittest
 import sys
+import os
 from SimPEG import Examples
 import numpy as np
 
+class compareInitFiles(unittest.TestCase):
+    def test_compareInitFiles(self):
+        print 'Checking that __init__.py up-to-date in SimPEG/Examples'
+        fName = os.path.abspath(__file__)
+        ExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['SimPEG', 'Examples'])
+
+        files = os.listdir(ExamplesDir)
+
+        pyfiles = []
+        [pyfiles.append(py.rstrip('.py')) for py in files if py.endswith('.py') and py != '__init__.py']
+
+        setdiff = set(pyfiles) - set(Examples.__examples__)
+
+        print ' Any missing files? ', setdiff 
+
+        didpass = (setdiff == set())
+
+        self.assertTrue(didpass, "Examples not up to date, run 'python __init__.py' from SimPEG/Examples to update")
+
 def get(test):
     def test_func(self):
         print '\nTesting %s.run(plotIt=False)\n'%test
@@ -10,11 +30,11 @@ def get(test):
         self.assertTrue(True)
     return test_func
 attrs = dict()
+
 for test in Examples.__examples__:
     attrs['test_'+test] = get(test)
 
 TestExamples = type('TestExamples', (unittest.TestCase,), attrs)
 
-
 if __name__ == '__main__':
     unittest.main()

tests/utils/test_Zero.py

@@ -1,5 +1,5 @@
 import unittest
-from SimPEG.Utils import Zero, Identity, sdiag
+from SimPEG.Utils import Zero, Identity, sdiag, mkvc
 from SimPEG import np, sp
 
 class Tests(unittest.TestCase):
@@ -29,6 +29,11 @@ class Tests(unittest.TestCase):
         assert a == 1
         self.assertRaises(ZeroDivisionError, lambda:3/z)
 
+        assert mkvc(z) == 0
+        assert sdiag(z)*a == 0
+        assert z.T == 0
+        assert z.transpose == 0
+
     def test_mat_zero(self):
         z = Zero()
         S = sdiag(np.r_[2,3])