diff --git a/SimPEG/MT/Problem1D/Probs.py b/SimPEG/MT/Problem1D/Probs.py
index 6232b2e9..418e78c2 100644
--- a/SimPEG/MT/Problem1D/Probs.py
+++ b/SimPEG/MT/Problem1D/Probs.py
@@ -23,7 +23,7 @@ class eForm_psField(BaseMTProblem):
     we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
 
     .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
+        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
     which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
 
     The primary field is estimated from a background model (commonly half space ).
@@ -40,6 +40,23 @@ class eForm_psField(BaseMTProblem):
         BaseMTProblem.__init__(self, mesh, **kwargs)
         self.fieldsPair = Fields1D_e
         # self._sigmaPrimary = sigmaPrimary
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
 
     @property
     def sigmaPrimary(self):
@@ -48,6 +65,7 @@ class eForm_psField(BaseMTProblem):
 
         """
         return self._sigmaPrimary
+
     @sigmaPrimary.setter
     def sigmaPrimary(self, val):
         # Note: TODO add logic for val, make sure it is the correct size.
@@ -62,16 +80,14 @@ class eForm_psField(BaseMTProblem):
             :return: A
         """
 
-        Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
         # Note: need to use the code above since in the 1D problem I want
         # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
         # Possible that _fieldType and _eqLocs can fix this
-        # Mmui = self.MfMui
-        # Msig = self.MeSigma
+        MeMui = self.MfMui
+        MfSigma = self.MfSigma
         C = self.mesh.nodalGrad
         # Make A
-        A = C.T*Mmui*C + 1j*omega(freq)*Msig
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
         # Either return full or only the inner part of A
         return A
 
@@ -81,15 +97,15 @@ class eForm_psField(BaseMTProblem):
         """
 
         dsig_dm = self.curModel.sigmaDeriv
-        MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        MeMui = self.MeMui
         #
         u_src = u['e_1dSolution']
-        dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
+        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
         if adjoint:
-            return 1j * omega(freq) * (  dMf_dsig.T * v )
+            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
         # Note: output has to be nN/nF, not nC/nE.
         # v should be nC
-        return 1j * omega(freq) * ( dMf_dsig * v )
+        return 1j * omega(freq) * ( dMfSigma_dm * v )
 
     def getRHS(self, freq):
         """
@@ -169,6 +185,23 @@ class eForm_TotalField(BaseMTProblem):
 
     def __init__(self, mesh, **kwargs):
         BaseMTProblem.__init__(self, mesh, **kwargs)
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
 
     def getA(self, freq, full=False):
         """
@@ -180,31 +213,24 @@ class eForm_TotalField(BaseMTProblem):
             :return: A
         """
 
-        Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
         # Note: need to use the code above since in the 1D problem I want
         # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
         # Possible that _fieldType and _eqLocs can fix this
-        # Mmui = self.MfMui
-        # Msig = self.MeSigma
+        # MeMui = self.MfMui
+        # MfSigma = self.MfSigma
         C = self.mesh.nodalGrad
         # Make A
-        A = C.T*Mmui*C + 1j*omega(freq)*Msig
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
         # Either return full or only the inner part of A
         if full:
             return A
         else:
             return A[1:-1,1:-1]
 
-    def getADeriv(self, freq, u, v, adjoint=False):
-        sig = self.curTModel
-        dsig_dm = self.curTModelDeriv
-        dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
-
-        if adjoint:
-            return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
-
-        return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        raise NotImplementedError('getADeriv is not implemented')
 
     def getRHS(self, freq):
         """
@@ -230,7 +256,7 @@ class eForm_TotalField(BaseMTProblem):
 
         return -Aio*eBC, eBC
 
-    def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
+    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
         raise NotImplementedError('getRHSDeriv not implemented yet')
         return None
 
diff --git a/SimPEG/MT/Utils/MT1Danalytic.py b/SimPEG/MT/Utils/MT1Danalytic.py
index cf0847d9..a28777cb 100644
--- a/SimPEG/MT/Utils/MT1Danalytic.py
+++ b/SimPEG/MT/Utils/MT1Danalytic.py
@@ -1,6 +1,7 @@
 # Analytic solution of EM fields due to a plane wave
 
 import numpy as np, SimPEG as simpeg
+from scipy.constants import mu_0, epsilon_0 as eps_0
 
 def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     '''Analytic solution for MT 1D layered earth. Returns E and H fields.
@@ -17,8 +18,8 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     # Note add an error check for the mesh and sigma are the same size.
 
     # Constants: Assume constant
-    mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
-    eps = 8.85*1e-12*np.ones((m1d.nC+1))
+    mu = mu_0*np.ones((m1d.nC+1))
+    eps = eps_0*np.ones((m1d.nC+1))
     # Angular freq
     w = 2*np.pi*freq
     # Add the halfspace value to the property
@@ -83,10 +84,6 @@ def getImpedance(m1d,sigma,freq):
 
     """
 
-    # Define constants
-    mu0 = 4*np.pi*1e-7
-    eps0 = 8.85e-12
-
     # Initiate the impedances
     Z1d = np.empty(len(freq) , dtype='complex')
     h = m1d.hx   #vectorNx[:-1]
@@ -95,13 +92,13 @@ def getImpedance(m1d,sigma,freq):
         om = 2*np.pi*fr
         Zall = np.empty(len(h)+1,dtype='complex')
         # Calculate the impedance for the bottom layer
-        Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
+        Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om)
 
         for nr,hi in enumerate(h):
             # Calculate the wave number
             # print nr,sigma[nr]
-            k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
-            Z = (mu0*om)/k
+            k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
+            Z = (mu_0*om)/k
 
             Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))