diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py
new file mode 100644
index 00000000..d01ae706
--- /dev/null
+++ b/SimPEG/EM/Static/SIP/ProblemSIP.py
@@ -0,0 +1,422 @@
+from SimPEG import Problem, Utils, Maps, Mesh
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from SimPEG.EM.Static.DC import getxBCyBC_CC
+from SurveySIP import Survey
+
+class ColeColePropMap(Maps.PropMap):
+    """
+        Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
+    """
+
+    eta = Maps.Property("Electrical Conductivity", defaultInvProp=True)
+    tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap))
+    taui   = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap))
+    c   = Maps.Property("Electrical Conductivity", defaultVal=1.)
+
+
+class BaseSIPProblem(BaseEMProblem):
+
+    surveyPair = Survey
+    fieldsPair = Fields
+    PropMap = ColeColePropMap
+    Ainv = None
+    sigma = None
+    rho = None
+    f = None
+    Ainv = None
+
+    def DebyeTime(t):
+        peta = self.curModel.eta*np.exp(-self.curModel.taui*t)
+        return peta
+
+    def EtaDeriv(t):
+        return np.exp(-self.curModel.taui*t)
+
+    def TauiDeriv(t):
+        return -self.curModel.eta*t*np.exp(-self.curModel.taui*t)
+
+    def fields(self, m):
+        self.curModel = m
+        if self.f is None:
+            self.f = self.fieldsPair(self.mesh, self.survey)
+            if self.Ainv == None:
+                A = self.getA()
+                self.Ainv = self.Solver(A, **self.solverOpts)
+            RHS = self.getRHS()
+            u = self.Ainv * RHS
+            Srcs = self.survey.srcList
+            self.f[Srcs, self._solutionType] = u
+        return self.f
+
+    def forward(self, m, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        Jv = self.dataPair(self.survey) #same size as the data
+        # A = self.getA()
+        JvAll = []
+        for tind in range(len(self.survey.times)):
+            #Pseudo-chareability
+            v = DebyeTime(self.survey.times[tind])
+            for src in self.survey.srcList:
+                u_src = f[src, self._solutionType] # solution vector
+                dA_dm_v = self.getADeriv(u_src, v)
+                dRHS_dm_v = self.getRHSDeriv(src, v)
+                du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+                for rx in src.rxList:
+                    timeindex = rx.getTimeP(self.survey.times)
+                    if timeindex[tind]:
+                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                        Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+            JvAll.append(Utils.mkvc(Jv))
+        # Conductivity (d u / d log sigma)
+        if self._formulation is 'EB':
+            return -np.hstack(JvAll)
+        # Conductivity (d u / d log rho)
+        if self._formulation is 'HJ':
+            return np.hstack(JvAll)
+
+    # def Jvec(self, m, v, f=None):
+
+    #     if f is None:
+    #         f = self.fields(m)
+
+    #     self.curModel = m
+
+    #     Jv = self.dataPair(self.survey) #same size as the data
+    #     x1 =
+    #     x2 =
+    #     # A = self.getA()
+    #     for src in self.survey.srcList:
+    #         u_src = f[src, self._solutionType] # solution vector
+
+    #         dA_dm_v = self.getADeriv(u_src, v)
+    #         dRHS_dm_v = self.getRHSDeriv(src, v)
+    #         du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+    #         for rx in src.rxList:
+
+    #             for tind in range(len(self.survey.times)):
+    #                 df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+    #                 df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+    #                 Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+    #     # Conductivity (d u / d log sigma)
+    #     if self._formulation is 'EB':
+    #         return -Utils.mkvc(Jv)
+    #     # Conductivity (d u / d log rho)
+    #     if self._formulation is 'HJ':
+    #         return Utils.mkvc(Jv)
+
+    # def Jtvec(self, m, v, f=None):
+    #     if f is None:
+    #         f = self.fields(m)
+
+    #     self.curModel = m
+
+    #     # Ensure v is a data object.
+    #     if not isinstance(v, self.dataPair):
+    #         v = self.dataPair(self.survey, v)
+
+    #     Jtv = np.zeros(m.size)
+    #     # AT = self.getA()
+
+    #     for src in self.survey.srcList:
+    #         u_src = f[src, self._solutionType]
+    #         for rx in src.rxList:
+    #             PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
+    #             df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+    #             df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
+    #             ATinvdf_duT = self.Ainv * df_duT
+    #             dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+    #             dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
+    #             du_dmT = -dA_dmT + dRHS_dmT
+    #             Jtv += df_dmT + du_dmT
+    #     # Conductivity ((d u / d log sigma).T)
+    #     if self._formulation is 'EB':
+    #         return -Utils.mkvc(Jtv)
+    #     # Conductivity ((d u / d log rho).T)
+    #     if self._formulation is 'HJ':
+    #         return Utils.mkvc(Jtv)
+
+    def getSourceTerm(self):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+    @property
+    def deleteTheseOnModelUpdate(self):
+        toDelete = []
+        return toDelete
+
+    # assume log rho or log cond
+    @property
+    def MeSigma(self):
+        """
+            Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
+        """
+        if getattr(self, '_MeSigma', None) is None:
+            self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
+        return self._MeSigma
+
+    @property
+    def MfRhoI(self):
+        """
+            Inverse of :code:`MfRho`
+        """
+        if getattr(self, '_MfRhoI', None) is None:
+            self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
+        return self._MfRhoI
+
+    def MfRhoIDeriv(self,u):
+        """
+            Derivative of :code:`MfRhoI` with respect to the model.
+        """
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
+        drho_dlogrho = Utils.sdiag(self.rho)
+        return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
+
+    # TODO: This should take a vector
+    def MeSigmaDeriv(self, u):
+        """
+            Derivative of MeSigma with respect to the model
+        """
+        dsigma_dlogsigma = Utils.sdiag(self.sigma)
+        return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
+
+class Problem3D_CC(BaseSIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseSIPProblem.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI D^\\top V
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        # TODO: this won't work for full anisotropy
+        MfRhoI = self.MfRhoI
+        A = D * MfRhoI * G
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * A
+        return A
+
+    def getADeriv(self, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        MfRhoIDeriv = self.MfRhoIDeriv
+
+        if adjoint:
+            # if self._makeASymmetric is True:
+            #     v = V * v
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
+        return D * (MfRhoIDeriv( G * u ) * v)
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return self.Vol.T * RHS
+
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+
+class Problem3D_N(BaseSIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # N potentials means B is on faces
+    fieldsPair    = Fields_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseSIPProblem.__init__(self, mesh, **kwargs)
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI D^\\top V
+
+        """
+
+        # TODO: this won't work for full anisotropy
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        A = Grad.T * MeSigma * Grad
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+
+        return A
+
+    def getADeriv(self, u, v, adjoint=False):
+        """
+
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        """
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        if not adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
+        elif adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
+
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+if __name__ == '__main__':
+
+
+    cs = 12.5
+    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hz = [(cs,7, -1.3),(cs,20)]
+    mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+    sigma = np.ones(mesh.nC)
+    prob = BaseSIPProblem(mesh, sigma=sigma)
+
+
diff --git a/SimPEG/EM/Static/SIP/RxSIP.py b/SimPEG/EM/Static/SIP/RxSIP.py
new file mode 100644
index 00000000..f4ff1a2f
--- /dev/null
+++ b/SimPEG/EM/Static/SIP/RxSIP.py
@@ -0,0 +1,82 @@
+import SimPEG
+import numpy as np
+from SimPEG.Utils import Zero, closestPoints
+
+class BaseRx(SimPEG.Survey.BaseTimeRx):
+    locs = None
+    rxType = None
+
+    knownRxTypes = {
+                    'phi':['phi',None],
+                    'ex':['e','x'],
+                    'ey':['e','y'],
+                    'ez':['e','z'],
+                    'jx':['j','x'],
+                    'jy':['j','y'],
+                    'jz':['j','z'],
+                    }
+
+    def __init__(self, locs, times, rxType, **kwargs):
+        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    def projGLoc(self, f):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        comp = self.knownRxTypes[self.rxType][1]
+        if comp is not None:
+            return f._GLoc(self.rxType) + comp
+        return f._GLoc(self.rxType)
+
+    def getTimeP(self, timesall):
+        """
+            Returns the time projection matrix.
+
+            .. note::
+
+                This is not stored in memory, but is created on demand.
+        """
+        time_inds = np.in1d(timesall, self.times)
+        return time_inds
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        if not adjoint:
+            return P*v
+        elif adjoint:
+            return P.T*v
+
+
+# DC.Rx.Dipole(locs)
+class Dipole(BaseRx):
+
+    def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        locs = [locsM, locsN]
+        # We may not need this ...
+        BaseRx.__init__(self, locs, times, rxType)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0] * len(self.times)
+
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+
+    def getP(self, mesh, Gloc):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
+        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
+        P = P0 - P1
+
+        if self.storeProjections:
+            self._Ps[mesh] = P
+
+        return P
diff --git a/SimPEG/EM/Static/SIP/SurveySIP.py b/SimPEG/EM/Static/SIP/SurveySIP.py
new file mode 100644
index 00000000..de375794
--- /dev/null
+++ b/SimPEG/EM/Static/SIP/SurveySIP.py
@@ -0,0 +1,32 @@
+import SimPEG
+from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG import np, sp, Survey
+from SimPEG.Utils import Zero, Identity
+from SimPEG.EM.Static.DC.SrcDC import BaseSrc
+from SimPEG.EM.Static.SIP.RxSIP import BaseRx
+
+class Survey(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+    times = None
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+        self.getUniqueTimes()
+
+    def getUniqueTimes(self):
+        time_rx = []
+        for src in self.srcList:
+            for rx in src.rxList:
+                time_rx.append(rx.times)
+        self.times = np.unique(np.hstack(time_rx))
+
+    def dpred(self, m, f=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pf(m)
+        """
+        return self.prob.Jvec(m, m, f=f)
diff --git a/SimPEG/EM/Static/SIP/__init__.py b/SimPEG/EM/Static/SIP/__init__.py
new file mode 100644
index 00000000..f81ceb3c
--- /dev/null
+++ b/SimPEG/EM/Static/SIP/__init__.py
@@ -0,0 +1,2 @@
+from ProblemSIP import Problem3D_CC, Problem3D_N
+from SurveySIP import Survey
diff --git a/SimPEG/EM/Static/__init__.py b/SimPEG/EM/Static/__init__.py
index c9b4dc0d..f6e0e757 100644
--- a/SimPEG/EM/Static/__init__.py
+++ b/SimPEG/EM/Static/__init__.py
@@ -1,2 +1,3 @@
 import DC
 import IP
+import SIP