diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index 46576df5..cab6c02f 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -239,14 +239,51 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 
 
 
-# class UpdateReferenceModel(Parameter):
+class update_IRLS(InversionDirective):
 
-#     mref0 = None
+     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.)
 
-#     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
+         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))
+#==============================================================================
diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py
new file mode 100644
index 00000000..c3517360
--- /dev/null
+++ b/SimPEG/Examples/DC_PseudoSection_Simulation.py
@@ -0,0 +1,179 @@
+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()
\ No newline at end of file
diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
index 00bf0259..c5754c93 100644
--- a/SimPEG/Mesh/DiffOperators.py
+++ b/SimPEG/Mesh/DiffOperators.py
@@ -565,7 +565,58 @@ 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
diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 4f2cb062..20a74d0d 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -989,5 +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) )           
+    
+                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
\ No newline at end of file
diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index f24b57de..f8d8257c 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -213,5 +213,20 @@ class BaseTimeProblem(BaseProblem):
         if hasattr(self, '_timeMesh'):
             del self._timeMesh
 
+class LinearProblem(BaseProblem):
+    
+    surveyPair = Survey.LinearSurvey
 
+    def __init__(self, mesh, G, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, **kwargs)
+        self.G = G
+
+    def fields(self, m):
+        return self.G.dot(m)
+
+    def Jvec(self, m, v, u=None):
+        return self.G.dot(v)
+
+    def Jtvec(self, m, v, u=None):
+        return self.G.T.dot(v)
 
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index 677fba5c..cdcfd59a 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -282,3 +282,239 @@ 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
diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py
index 47a88ae2..60ef4330 100644
--- a/SimPEG/Survey.py
+++ b/SimPEG/Survey.py
@@ -1,6 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
 
-
 class BaseRx(object):
     """SimPEG Receiver Object"""
 
@@ -375,3 +374,11 @@ 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]
diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py
index 52d13820..435856f7 100644
--- a/SimPEG/Utils/ModelBuilder.py
+++ b/SimPEG/Utils/ModelBuilder.py
@@ -118,6 +118,44 @@ 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.