remove import of ipywidgets

# Conflicts:
#	SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
#	SimPEG/Examples/__init__.py
#	SimPEG/Examples/sphereElectrostatic_example.py
#	SimPEG/Optimization.py
#	docs/examples/DC_PseudoSection_Simulation.rst
#	docs/examples/Inversion_IRLS.rst
#	docs/examples/MT_1D_analytic_nlayer_Earth.rst

SimPEG/DCIP/DCIPUtils.py

@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):
 
     return model
 
-def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
     from scipy.interpolate import griddata
     import pylab as plt
 
-    # Set depth to 0 for now
-    z0 = 0.
-
     # Pre-allocate
     midx = []
     midz = []
@@ -259,38 +256,53 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
         midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
 
-    ax = axs
-
     # Grid points
     grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
     grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
-
+    
+    # Scale the color scheme
     if clim == None:
         vmin, vmax = rho.min(), rho.max()
     else:
         vmin, vmax = clim[0], clim[1]
 
+    # Plot data
     grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
-    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax))
-    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
-
-    cmin,cmax = cbar.get_clim()
-    ticks = np.linspace(cmin,cmax,3)
-    cbar.set_ticks(ticks)
-    cbar.ax.tick_params(labelsize=10)
+    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax)
+    plt.gca().tick_params(axis='both', which='major', labelsize=8)
     
-    if dtype == 'appc':
-        cbar.set_label("App.Cond",size=12)
-    elif dtype == 'appr':
-        cbar.set_label("App.Res.",size=12)
-    elif dtype == 'volt':
-        cbar.set_label("Potential (V)",size=12)
+    if contour is not None:
+        plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax)
+        
+    # Add scatter points
+    axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax)
+    
+    if colorbar:
+        
+        if dtype == 'volt':
+            cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
 
-    # Plot apparent resistivity
-    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
+        else:               
+            cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
+        
+        cmin,cmax = cbar.get_clim()
+        ticks = np.linspace(cmin,cmax,3)
+        cbar.set_ticks(ticks)
+        cbar.ax.tick_params(labelsize=10)
+        
+        if cblabel:
+            if dtype == 'appc':
+                cbar.set_label("App.Cond",size=12)
+            elif dtype == 'appr':
+                cbar.set_label("App.Res.",size=12)
+            elif dtype == 'volt':
+                cbar.set_label("Potential (V)",size=12)
 
-    #ax.set_xticklabels([])
-    #ax.set_yticklabels([])
+
+
+    if not axlabel:
+        axs.set_xticklabels([])
+        axs.set_yticklabels([])
 
     plt.gca().set_aspect('equal', adjustable='box')
 
@@ -448,15 +460,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
     """
         Write UBC GIF DCIP 2D or 3D observation file
 
         Input:
-        :string fileName -> including path where the file is written out
-        :DCsurvey -> DC survey class object
-        :string dtype ->  either '2D' | '3D'
-        :string  stype ->  either 'SURFACE' | 'GENERAL'
+        :string     fileName -> including path where the file is written out
+        :DCsurvey   DC survey class object
+        :string     dtype ->  either '2D' | '3D'
+        :string     stype ->  either 'SURFACE' | 'GENERAL'
 
         Output:
         :param UBC2D-Data file
@@ -471,10 +483,16 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
     assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
     assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-
+    
     fid = open(fileName,'w')
-    fid.write('! ' + stype + ' FORMAT\n')
-
+    
+    
+    if iptype!=0:
+        fid.write('IPTYPE=%i\n'%iptype)
+        
+    else:
+        fid.write('! ' + stype + ' FORMAT\n')
+        
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -498,7 +516,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
                 B = np.repeat(tx[0,1],M.shape[0],axis=0)
                 M = M[:,0]
                 N = N[:,0]
-
+                
                 np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
 
@@ -506,18 +524,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
                 if stype == 'SURFACE':
 
-                    fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
                     M = M[:,0]
                     N = N[:,0]
 
                 if stype == 'GENERAL':
 
+                    # Flip sign for z-elevation to depth
+                    tx[2::2,:] = -tx[2::2,:]
+                    
                     fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
                     M = M[:,0::2]
                     N = N[:,0::2]
 
+                    # Flip sign for z-elevation to depth
+                    M[:,1::2] = -M[:,1::2]
+                    N[:,1::2] = -N[:,1::2]
+                    
                 fid.write('%i\n'% nD)
-                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
 
         if dtype=='3D':
 
@@ -529,11 +554,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
             if stype == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx))
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
 
             fid.write('%i\n'% nD)
             np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
-
+            fid.write('\n')
+            
         count += nD
 
     fid.close()
@@ -640,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     DCsurvey2D.std = np.asarray(DCsurvey.std)
 
     return DCsurvey2D
-
-def readUBC_DC3Dobs(fileName):
+    
+def readUBC_DC3Dobs(fileName, dtype = 'DC'):
     """
-        Read UBC GIF DCIP 3D observation file and generate survey
+        Read UBC GIF IP 3D observation file and generate survey
 
         Input:
         :param fileName, path to the UBC GIF 3D obs file
 
         Output:
-        :param DCIPsurvey
+        :param IPsurvey
         :return
 
-        Created on Mon April 6th, 2015
-
         @author: dominiquef
 
     """
-
+    zflag = True # Flag for z value provided
+    
     # Load file
-    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-
+    if dtype == 'IP':
+        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
+        
+    elif dtype == 'DC':
+        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+        
+    else:
+        print "dtype must be 'DC'(default) | 'IP'"
+        
     # Pre-allocate
     srcLists = []
     Rx = []
     d = []
     wd = []
-    zflag = True # Flag for z value provided
+    
 
     # Countdown for number of obs/tx
     count = 0
     for ii in range(obsfile.shape[0]):
 
+        # Skip if blank line
         if not obsfile[ii]:
             continue
 
-        # First line is transmitter with number of receivers
+        # First line or end of a transmitter block, read transmitter info
         if count==0:
-
-            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+            # Read the line
+            temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
             count = int(temp[-1])
 
             # Check if z value is provided, if False -> nan
             if len(temp)==5:
-                tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
-                zflag = False
+                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
+                
+                zflag = False # Pass on the flag to the receiver loc
 
             else:
                 tx = temp[:-1]
@@ -692,8 +726,16 @@ def readUBC_DC3Dobs(fileName):
             rx = []
             continue
 
-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
 
+        # Filter out negative IP
+#        if temp[-2] < 0:       
+#            count = count -1
+#            print "Negative!"
+#            
+#        else:
+            
+        # If the Z-location is provided, otherwise put nan
         if zflag:
 
             rx.append(temp[:-2])
@@ -703,7 +745,7 @@ def readUBC_DC3Dobs(fileName):
                 wd.append(temp[-1])
 
         else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
             # Check if there is data with the location
             if len(temp)==6:
                 d.append(temp[-2])
@@ -711,7 +753,7 @@ def readUBC_DC3Dobs(fileName):
 
         count = count -1
 
-        # Reach the end of transmitter block
+        # Reach the end of transmitter block, append the src, rx and continue
         if count == 0:
             rx = np.asarray(rx)
             Rx = DC.RxDipole(rx[:,:3],rx[:,3:])

SimPEG/Directives.py

@@ -222,7 +222,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
             mref = 0
         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
         phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim==2:
+        if self.prob.mesh.dim >= 2:
             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
             phi_my = 0.5 * my.dot(my)
         else:
@@ -237,41 +237,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
-#==============================================================================
-# class SaveOutputDictEveryIteration(_SaveEveryIteration):
-#     """SaveOutputDictEveryIteration
-#     A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
-#     """
-#
-#     def initialize(self):
-#         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
-#
-#     def endIter(self):
-#         # Save the data.
-#         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-#         phi_ms = 0.5*ms.dot(ms)
-#         if self.reg.mrefInSmooth == True:
-#             mref = self.reg.mref
-#         else:
-#             mref = 0
-#         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#         phi_mx = 0.5 * mx.dot(mx)
-#         if self.prob.mesh.dim==2:
-#             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_my = 0.5 * my.dot(my)
-#         else:
-#             phi_my = 'NaN'
-#         if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
-#             mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_mz = 0.5 * mz.dot(mz)
-#         else:
-#             phi_mz = 'NaN'
-#
-#
-#         # Save the file as a npz
-#         np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
-#
-#==============================================================================
 
 # class UpdateReferenceModel(Parameter):
 
@@ -293,6 +258,7 @@ class Update_IRLS(InversionDirective):
     phi_m_last = None
     phi_d_last = None
 
+
     def initialize(self):
 
         # Scale the regularization for changes in norm
@@ -310,7 +276,7 @@ class Update_IRLS(InversionDirective):
             self.phi_d_last = self.invProb.phi_d
 
     def endIter(self):
-        # Cool the threshold parameter
+        # Cool the threshold parameter if required
         if getattr(self, 'factor', None) is not None:
             eps = self.reg.eps / self.factor
 
@@ -325,28 +291,44 @@ class Update_IRLS(InversionDirective):
          # Update the model used for the IRLS weights
         self.reg.curModel = self.invProb.curModel
 
-        # Temporarely set gamma to 1.
+        # Temporarely set gamma to 1. to get raw phi_m
         self.reg.gamma = 1.
 
-        # Compute change in model objective function and update scaling
+        # Compute new model objective function value
         phim_new = self.reg.eval(self.invProb.curModel)
 
+        # Update gamma to scale the regularization between IRLS iterations
         self.reg.gamma = self.phi_m_last / phim_new
 
-        self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+        # Set the weighting matrix to None so that it is recomputed next time
+        # it is called in the inversion
+        self.reg._W = None
 
 class Update_lin_PreCond(InversionDirective):
-
-
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
+    
+    def initialize(self):
+    
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # 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.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+            
     def endIter(self):
         # Cool the threshold parameter
-
+        if self.onlyOnStart==True:
+            return
+            
         if getattr(self.opt, 'approxHinv', None) is not None:
             # 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.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
+            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.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
-            print 'Updated pre-cond'
+
 
 class Update_Wj(InversionDirective):
     """
@@ -373,3 +355,19 @@ class Update_Wj(InversionDirective):
             JtJdiag = JtJdiag / max(JtJdiag)
 
             self.reg.wght = JtJdiag
+
+class Scale_Beta(InversionDirective):
+    """
+        Instead of a linear cooling schedule, beta is allowed to change based
+        on the ratio between the target misfit and the current data misfit. The
+        update is done only if the misfit is outside some threshold bounds.
+    """
+    tol = 0.05
+
+    def endIter(self):
+
+        # Check if misfit is within the tolerance, otherwise adjust beta
+        val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+        if np.abs(1.-val) > self.tol:
+            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

SimPEG/Examples/Inversion_IRLS.py

@@ -86,12 +86,12 @@ def run(N=200, plotIt=True):
     #reg.recModel = mrec
     reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 2e-3
-    reg.eps_q = 2e-3
+    reg.eps_p = 5e-2
+    reg.eps_q = 1e-2
     reg.norms   = [0., 0., 2., 2.]
     reg.wght = wr
 
-    opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3)
+    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
     invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
     beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
     #betaest = Directives.BetaEstimate_ByEig()

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -0,0 +1,427 @@
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+import numpy as np
+from SimPEG.EM.Utils import k, omega
+
+"""
+MT1D: n layered earth problem
+*****************************
+
+Author: Thibaut Astic
+Contact: thast@eos.ubc.ca
+
+
+This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
+
+We start by looking at Maxwell's equations in the electric
+field \\\(\\\mathbf{E}\\) and the magnetic flux
+\\\(\\\mathbf{H}\\) to write the wave equations
+\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
+\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
+
+Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
+\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
+\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
+
+With U and D the Up and Down components of the E-field.
+
+The iteration from one layer to another is ensure by:
+
+\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
+    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
+
+And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
+and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
+
+The layer 0 is assumed to be the air layer.
+
+"""
+
+#Define a frquency range for a survey
+frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
+
+#Functions to create random physical Properties for a n-layered earth
+thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
+                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
+                                                            ,decimals =1))
+
+sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
+                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
+                                                      ,decimals=3))
+
+mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
+                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
+                                                    ,decimals=1))
+
+eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
+                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
+                                                                ,decimals=1))
+
+#Evaluate Impedance Z of a layer
+ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
+
+#Complex Cole-Cole Conductivity - EM utils
+PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
+
+#Converted thickness array into top of layer array
+top = lambda thick: np.cumsum(thick)
+
+#Propagation Matrix and theirs inverses
+
+#matrix T for transition of Up and Down components accross a layer
+T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
+
+Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
+
+#transition of Up and Down components accross a layer
+UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
+
+
+#matrix P relating Up and Down components with E and H fields
+P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
+
+Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
+
+
+#Time Variation of E and H
+E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
+H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
+
+#Plot the configuration of the problem
+def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
+
+    topn = top(thick)
+    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
+
+    ax.set_ylim([z.min(),z.max()])
+    ax.set_xlim([-widthg,widthg])
+
+    ax.set_ylabel("Depth (m)", fontsize=16.)
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+    #define filling for the different layers
+    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
+
+    #Write the physical properties of air
+    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
+            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
+    for i in range(1,len(topn)-1,1):
+        if topn[i] == topn[i+1]:
+            pass
+        else:
+            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
+                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
+                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
+                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
+            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
+
+    #Write the physical properties of the n-th layer and fill it with pattern
+    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
+    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
+
+    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
+            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
+            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
+            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    #plot Trees!
+    ax.annotate("",
+            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+
+    ax.invert_yaxis()
+
+    return ax
+
+#Propagate Up and Down component for a certain frequency & evaluate E and H field
+
+def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
+
+    sigcm = np.zeros_like(sig,dtype='complex_')
+
+    for j in range(1,len(sig)):
+        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
+
+    K = k(f, sigcm, mu, eps)
+    Z = ImpZ(f,mu,K)
+
+    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+
+    UD[1,-1] = 1.
+
+    for i in range(-2,-(n+2),-1):
+
+        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
+        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
+
+    for j in range(0,n+1):
+        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
+
+    return UD, EH, Z ,K
+
+
+#Evaluate the apparent resistivity and phase for a frequency range
+def appres(F,H,sig,chg,taux,c,mu,eps,n):
+
+    Res = np.zeros_like(F)
+    Phase = np.zeros_like(F)
+    App_ImpZ= np.zeros_like(F,dtype='complex_')
+
+    for i in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
+
+        App_ImpZ[i] = EH[0,1]/EH[1,1]
+
+        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
+        Phase[i] = np.angle(App_ImpZ[i], deg = True)
+
+    return Res,Phase
+
+#Evaluate Up, Down components, E and H field, for a frequency range,
+#a discretized depth range and a time range (use to calculate envelope)
+def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
+
+    topc = top(H)
+
+    layer = np.zeros(len(zsample),dtype=np.int)-1
+
+    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
+
+    for i in range(0,n+1,1):
+        layer = layer+(zsample>=topc[i])*1
+
+    for j in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
+
+        for p in range(0,len(zsample)):
+
+            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
+
+            for q in range(0,len(tsample)):
+
+                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
+                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
+                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+
+    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
+
+
+#Function to Plot Apparent Resistivity and Phase
+def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
+
+    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
+
+    fig,ax = plt.subplots(1,2,figsize=(16,10))
+
+    ax[0].scatter(Res,F,color='black')
+    ax[0].set_xscale('Log')
+    ax[0].set_yscale('Log')
+    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
+    ax[0].set_ylim([F.min(),F.max()])
+    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
+    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
+    ax[0].grid(which='major')
+
+    ax0 = ax[0].twiny()
+
+    ax0.set_xlim([0.,90.])
+    ax0.set_ylim([F.min(),F.max()])
+    ax0.scatter(Phase,F,color='purple')
+    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
+
+    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
+
+    ax[0].tick_params(labelsize=16)
+    ax[1].tick_params(labelsize=16)
+    ax0.tick_params(labelsize=16)
+
+    if PlotEnvelope:
+
+        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
+        fenvelope1n=np.ones(100)*fenvelope
+        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
+
+        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
+        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
+
+        ax1=ax[1].twiny()
+
+        ax[1].tick_params(labelsize=16)
+        ax1.tick_params(labelsize=16)
+
+        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
+
+        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
+
+        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
+                      color='blue', alpha=0.1)
+
+        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
+                      color='red', alpha=0.1)
+
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+    else:
+        print 'No envelop (if True, might be slow)'
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
+        ax[1].get_xaxis().set_ticks([])
+
+    plt.show()
+
+#Interactive MT for Notebook
+def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
+
+    frangn=frange(-5,5,100.)
+    sig3= np.array([0.,0.001,0.1, 0.001])
+    thick3 = np.array([120000.,50.,50.])
+    eps3=np.array([1.,1.,1.,1])
+    mu3=np.array([1.,1.,1.,1])
+    chg3=np.array([0.,0.1,0.,0.2])
+    chg3_0=np.array([0.,0.1,0.,0.])
+    taux3=np.array([0.,0.1,0.,0.1])
+    c3=np.array([1.,1.,1.,1.])
+
+    sig3[1]=sigl1
+    sig3[1]=10.**sig3[1]
+    sig3[2]=sigl2
+    sig3[2]=10.**sig3[2]
+    sig3[3]=sigl3
+    sig3[3]=10.**sig3[3]
+    mu3[1]=mul1
+    mu3[2]=mul2
+    mu3[3]=mul3
+    eps3[1]=epsl1
+    eps3[2]=epsl2
+    eps3[3]=epsl3
+    thick3[1]=h1
+    thick3[2]=h2
+
+    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
+
+
+def run(plotIt=True, n=3):
+    # something to make a plot
+
+    F = frange(-5.,5.,20)
+    H = thick(50.,100.,n)
+    sign = sig(-5.,0.,n)
+    mun = mu(1.,2.,n)
+    epsn = eps(1.,9.,n)
+    chg = np.zeros_like(sign)
+    taux = np.zeros_like(sign)
+    c = np.zeros_like(sign)
+
+    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
+
+    if plotIt:
+
+        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
+
+    return Res, Phase
+
+if __name__ == '__main__':
+    run(plotIt=True)
+
+
+
+
+

SimPEG/Examples/__init__.py

@@ -18,10 +18,12 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
+import MT_1D_analytic_nlayer_Earth
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
+import sphereElectrostatic_example
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "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_ForwardAndInversion", "MT_3D_Foward"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "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", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "sphereElectrostatic_example"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Examples/sphereElectrostatic_example.py

@@ -0,0 +1,785 @@
+from scipy.constants import epsilon_0
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+import numpy as np
+from SimPEG.Utils import ndgrid, mkvc
+
+'''
+Authors: Thibaut Astic, Lindsey Heagy, Sanna Tyrvainen, Ronghua Peng
+
+
+This code defines function to resolve analytically the electrostatic sphere problem.
+We first define a problem configuration, with a conductive or resistive sphere in a
+wholespace background.
+We then calculate the potential, then the electric field, then the current density and
+finally the charges accumulation.
+
+Several plotting functions are defined for data visualisation.
+
+
+'''
+
+# Plot options
+ftsize_title = 18      #font size for titles
+ftsize_axis  = 14      #font size for axis ticks
+ftsize_label = 14      #font size for axis labels
+
+# Radius function, useful sigma ratio, and log scale converter
+r  = lambda x,y,z: np.sqrt(x**2.+y**2.+z**2.)
+sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0)
+
+#tools to convert log conductivity in conductivity
+def conductivity_log_wrapper(log_sig0,log_sig1):
+    sig0 = 10.**log_sig0
+    sig1 = 10.**log_sig1
+
+    return sig0,sig1
+
+# Examples
+#Plot the configuration. Label=False is used to generate a general case figure
+def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
+    '''
+    XYZ: ndgrid
+    sig0: conductivity of the background
+    sig1: conductivity of the sphere
+    R: radius of the sphere
+    E0: Amplitude of the uniform electrostatic field
+    ax: ax where to plot the configuration
+    label: True: plot real values, False: plot general case
+    colorsphere: color of the sphere, format [x,x,x]
+    '''
+
+    xplt = np.linspace(-R, R, num=100)
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    dx = xr[1]-xr[0]
+    top = np.sqrt(R**2-xplt**2)
+    bot = -np.sqrt(R**2-xplt**2)
+
+    if R != 0:
+        ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5)
+        ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 )
+        ax.arrow(0.,0.,np.sqrt(2.)*R/2.,np.sqrt(2.)*R/2.,head_width=0.,head_length=0.)
+
+        if label:
+            ax.annotate(("$\sigma_1$=%3.3f mS/m")%(sig1*10.**(3.)),
+            xy=(0.,-R/2.), xycoords='data',
+            xytext=(0.,-R/2.), textcoords='data',
+            fontsize=14.)
+            ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
+            xy=(0.,-1.5*R), xycoords='data',
+            xytext=(0.,-1.5*R), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$R$ = %1i m')%(R),
+            xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
+            xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
+            fontsize=14.)
+            ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+            ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+            ax.tick_params(labelsize=ftsize_axis)
+
+        else:
+            ax.set_xticklabels([])
+            ax.set_yticklabels([])
+            ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14)
+            ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14)
+            ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$R$'),
+            xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
+            xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
+            fontsize=14.)
+            ax.set_xlabel('x',fontsize=12)
+            ax.set_ylabel('y',fontsize=12)
+
+    else:
+        if label:
+            ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
+            xy=(0.,-1.5*R), xycoords='data',
+            xytext=(0.,-1.5*R), textcoords='data',
+            fontsize=14.)
+            ax.annotate(('$\mathbf{E_0} = %1i  \mathbf{\hat{x}}$ V/m')%(E0),
+            xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
+            xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
+            fontsize=14.)
+            ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+            ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+            ax.tick_params(labelsize=ftsize_axis)
+
+        else:
+            ax.set_xticklabels([])
+            ax.set_yticklabels([])
+            ax.text(-0.05,-10,'$\sigma_0$',fontsize=14)
+            ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14)
+            ax.set_xlabel('x',fontsize=12)
+            ax.set_ylabel('y',fontsize=12)
+
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    [ax.arrow(xr.min(),_,np.abs(xr.max()-xr.min())/20.,0.,head_width=5.,head_length=2.,color='k') for _ in np.linspace(yr.min(),yr.max(),num=10)]
+    ax.patch.set_facecolor([0.4,0.7,0.4])
+    ax.patch.set_alpha(0.2)
+
+    ax.set_aspect('equal')
+
+
+
+    return ax
+
+def get_Conductivity(XYZ,sig0,sig1,R):
+    '''
+    Define the conductivity for each point of the space
+    '''
+    x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
+    r_view=r(x,y,z)
+
+    ind0= (r_view>R)
+    ind1= (r_view<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Sigma = np.zeros_like(x)
+
+    Sigma[ind0] = sig0
+    Sigma[ind1] = sig1
+
+    return Sigma
+
+
+def get_Potential(XYZ,sig0,sig1,R,E0):
+
+    '''
+    Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin
+    :input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
+    '''
+
+    x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
+
+    sig_cur = sigf(sig0,sig1)
+
+    r_cur = r(x,y,z)  # current radius
+
+    ind0 = (r_cur > R)
+    ind1 = (r_cur <= R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Vt = np.zeros_like(x)
+    Vp = np.zeros_like(x)
+    Vs = np.zeros_like(x)
+
+    Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere
+    Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0)             # inside the sphere
+
+
+    Vp = - E0*x  # primary potential
+
+    Vs = Vt - Vp # secondary potential
+
+    return Vt,Vp,Vs
+
+#plot the primary potential on ax
+def Plot_Primary_Potential(XYZ,Vp,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Primary Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+#plot the total potential on ax
+def Plot_Total_Potential(XYZ,Vt,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Total Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+#plot the secondary potential on ax
+def Plot_Secondary_Potential(XYZ,Vs,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size))
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Secondary Potential',fontsize=ftsize_title)
+    cb = plt.colorbar(Pplot,ax=ax)
+    cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_aspect('equal')
+    ax.tick_params(labelsize=ftsize_axis)
+
+    return ax
+
+
+def get_ElectricField(XYZ,sig0,sig1,R,E0):
+    '''
+    Function that returns the total, the primary and the secondary electric fields,
+    input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    r_cur=r(x,y,z)  # current radius
+
+    ind0= (r_cur>R)
+    ind1= (r_cur<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Ep = np.zeros(shape=(len(x),3))
+    Ep[:,0] = E0
+
+    Et = np.zeros(shape=(len(x),3))
+
+    Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.);
+    Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1);
+    Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1);
+
+    Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0;
+    Et[ind1,1] = 0.;
+    Et[ind1,2] = 0.;
+
+    Es = Et - Ep
+
+    return Et, Ep, Es
+
+#plot the total electric field on ax
+def Plot_Total_ElectricField(XYZ,Et,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    EtXr = Et[:,0].reshape(xr.size, yr.size)
+    EtYr = Et[:,1].reshape(xr.size, yr.size)
+    EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Eplot = ax.pcolor(xr,yr,EtAmp)
+    cb = plt.colorbar(Eplot,ax=ax)
+    cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05)
+    ax.set_title('Total Field',fontsize=ftsize_title)
+
+
+    return ax
+
+#plot the secondary electric field on ax
+def Plot_Secondary_ElectricField(XYZ,Es,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    xcirc = xr[np.abs(xr) <= R]
+
+    EsXr = Es[:,0].reshape(xr.size, yr.size)
+    EsYr = Es[:,1].reshape(xr.size, yr.size)
+    EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Eplot = ax.pcolor(xr,yr,EsAmp)
+    cb = plt.colorbar(Eplot,ax=ax)
+    cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_title('Secondary Field',fontsize=ftsize_title)
+
+    return ax
+
+
+def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es):
+    '''
+    Function that returns the total, the primary and the secondary current densities,
+    :input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields,
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    r_cur=r(x,y,z)
+
+    ind0= (r_cur>R)
+    ind1= (r_cur<=R)
+
+    assert (ind0 + ind1).all(), 'Some indicies not included'
+
+    Jt = np.zeros(shape=(len(x),3))
+    J0 = np.zeros(shape=(len(x),3))
+    Js = np.zeros(shape=(len(x),3))
+
+
+    Jp = sig0*Ep
+
+    Jt[ind0,:] = sig0*Et[ind0,:]
+    Jt[ind1,:] = sig1*Et[ind1,:]
+
+    Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:])
+    Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:]
+
+    return Jt,Jp,Js
+
+#plot the total currents density on ax
+def Plot_Total_Currents(XYZ,Jt,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    JtXr = Jt[:,0].reshape(xr.size, yr.size)
+    JtYr = Jt[:,1].reshape(xr.size, yr.size)
+    JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size))
+    cb = plt.colorbar(Jplot,ax=ax)
+    cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
+    ax.set_title('Total Current Density',fontsize=ftsize_title)
+
+    return ax
+
+
+#plot the secondary currents density on ax
+def Plot_Secondary_Currents(XYZ,Js,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    JsXr = Js[:,0].reshape(xr.size, yr.size)
+    JsYr = Js[:,1].reshape(xr.size, yr.size)
+    JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Js[:,2]**2).reshape(xr.size,yr.size)
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_aspect('equal')
+
+    Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size))
+    cb = plt.colorbar(Jplot,ax=ax)
+    cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
+    cb.ax.tick_params(labelsize=ftsize_axis)
+    ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
+    ax.set_title('Secondary Current Density',fontsize=ftsize_title)
+
+    return ax
+
+
+def get_ChargesDensity(XYZ,sig0,sig1,R,Ep):
+    '''
+    Function that returns the charges accumulation at the background/sphere interface,
+    :input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields,
+    '''
+
+    x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
+
+    dx = x[1]-x[0]
+
+    r_cur=r(x,y,z)
+
+    ind0 = (r_cur > R)
+    ind1 = (r_cur < R)
+    ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) )
+
+    assert (ind0 + ind1 + ind2).all(), 'Some indicies not included'
+
+    rho = np.zeros_like(x)
+
+    rho[ind0] = 0
+    rho[ind1] = 0
+    rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.))
+
+    return rho
+
+#Plot charges density on ax
+def Plot_ChargesDensity(XYZ,rho,R,ax):
+
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+    xcirc = xr[np.abs(xr) <= R]
+
+    ax.set_xlim([xr.min(),xr.max()])
+    ax.set_ylim([yr.min(),yr.max()])
+    ax.set_aspect('equal')
+    Cplot = ax.pcolor(xr,yr,rho.reshape(xr.size, yr.size))
+    cb1 = plt.colorbar(Cplot,ax=ax)
+    cb1.set_label(label= 'Charge Density ($C/m^2$)',size=ftsize_label) #weight='bold')
+    cb1.ax.tick_params(labelsize=ftsize_axis)
+    ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
+    ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
+    ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
+    ax.tick_params(labelsize=ftsize_axis)
+    ax.set_title('Charges Density', fontsize=ftsize_title)
+
+    return ax
+
+def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
+
+    '''
+    Function that return array of midpoints electrodes, electrodes positions,
+    potentials differences for total and secondary potentials fields, unormalized and
+    normalized to electrodes distances.
+    sig0: background conductivity
+    sig1: sphere conductivity
+    R: Sphere's radius
+    E0: uniform E field value
+    start: start point for the profile start.shape = (2,)
+    end: end point for the profile end.shape = (2,)
+    nbmp: number of dipoles
+    mn: Space between the M and N electrodes
+    '''
+
+    #D: total distance from start to end
+    D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.)
+
+    #MP: dipoles'midpoint positions (x,y)
+    MP = np.zeros(shape=(nbmp,2))
+    MP[:,0] = np.linspace(start[0],end[0],nbmp)
+    MP[:,1] = np.linspace(start[1],end[1],nbmp)
+
+    #Dipoles'Electrodes positions around each midpoints
+    EL = np.zeros(shape=(2*nbmp,2))
+    for n in range(0,len(EL),2):
+        EL[n,0]   = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2.
+        EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2.
+        EL[n,1]   = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2.
+        EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2.
+
+    VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL)
+    VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL)
+    dVtMP = np.zeros(nbmp)  #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP)
+    dVtMPn = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+    dVsMP = np.zeros(nbmp)  #Diffence (d-) of Secondaty Potential (Vt-) at each dipole (-MP)
+    dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+    dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP)
+    dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
+
+    #Computing VtEL
+    for m in range(0,2*nbmp):
+        if (r(EL[m,0],EL[m,1],0) > R):
+            VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.)
+        else:
+            VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0)
+
+    #Computing VsEL
+    VsEL = VtEL + E0*EL[:,0]
+
+    #Computing dVtMP, dVsMP
+    for p in range(0,nbmp):
+        dVtMP[p] = VtEL[2*p]-VtEL[2*p+1]
+        dVtMPn[p] = dVtMP[p]/mn
+        dVsMP[p] = VsEL[2*p]-VsEL[2*p+1]
+        dVsMPn[p] = dVsMP[p]/mn
+
+    return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn
+
+#Compare the DC response of two configurations
+def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor):
+
+    #Define the mesh
+    xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
+
+    #Defining the Profile
+    start = np.array([xstart,ystart])
+    end = np.array([xend,yend])
+
+    #Calculating the data from the defined survey line for Configuration 0 and 1
+    MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing)
+    MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing)
+
+
+    # Initializing the figure
+    fig = plt.figure(figsize=(20,20))
+    ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6)
+    ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6)
+    ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4)
+    ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6)
+    ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6)
+
+    #Plotting the Configuration 0
+    ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1])
+
+    #Plotting the Configuration 1
+    ax1   = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6])
+
+    #Plotting the Data (Legends)
+    ax2.set_title('Potential Differences',fontsize=ftsize_title)
+    ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
+    ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
+    ax2.tick_params(labelsize=ftsize_axis)
+    ax2.grid()
+
+    #Calculating the potential
+    Vt0,Vp0,Vs0 = get_Potential(XYZ,sig0,sig1,R0,E0)
+    Vt1,Vp1,Vs1 = get_Potential(XYZ,sig0,sig2,R1,E0)
+
+    if PlotOpt == 'Total':
+        ax3= Plot_Total_Potential(XYZ,Vt0,R0,ax3)
+        ax4= Plot_Total_Potential(XYZ,Vt1,R1,ax4)
+
+        #Plot the Data (from Configuration 0)
+        gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0
+                         ,marker='o',color='blue',linewidth=3.,label ='Left Model Response' )
+
+        #Plot the Data (from Configuration 1)
+        gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1
+                ,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
+        ax2.legend(('Left Model Response','Right Model Response'),loc=4)
+
+    elif PlotOpt == 'Secondary':
+        #plot the secondary potentials
+        ax3= Plot_Secondary_Potential(XYZ,Vt0,R0,ax3)
+        ax4= Plot_Secondary_Potential(XYZ,Vt1,R1,ax3)
+
+        #Plot the data(from configuration 0)
+        gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue'
+                ,marker='o',linewidth=3.,label ='Left Model Response' )
+
+
+        #Plot the Data (from Configuration 1)
+        gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1
+                 ,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
+        ax2.legend(('Left Model Response','Right Model Response'),loc=4 )
+
+    else:
+        print('What dont you get? Total or Secondary?')
+
+    #Legends
+    ax3.plot(MP0[:,0],MP0[:,1],color='gray')
+    Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black')
+    Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red')
+    ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
+
+    ax4.plot(MP1[:,0],MP1[:,1],color='gray')
+    Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black')
+    Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red')
+    ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
+
+    return fig
+
+#Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field
+def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b):
+
+    sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
+    E0   = 1.           # inducing field strength in V/m
+    n = 100             #level of discretisation
+    xr = np.linspace(-200., 200., n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+
+    Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
+
+
+    fig, ax = plt.subplots(1,2,figsize=(18,6))
+
+    #Setup figure 1 with options Configuration, Total or Secondary,
+    #then Potential, ElectricField, Current Density or Charges Density
+    if Figure1a == 'Configuration':
+        ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6])
+
+    elif Figure1a == 'Total':
+
+        if Figure1b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Total_Potential(XYZ,Vt,R,ax[0])
+
+        elif Figure1b == 'ElectricField':
+            ax[0] = Plot_Total_ElectricField(XYZ,Et,R,ax[0])
+
+        elif Figure1b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Total_Currents(XYZ,Jt,R,ax[0])
+
+        elif Figure1b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[0])
+
+    elif Figure1a == 'Secondary':
+
+        if Figure1b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Secondary_Potential(XYZ,Vs,R,ax[0])
+
+        elif Figure1b == 'ElectricField':
+            ax[0] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[0])
+
+        elif Figure1b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Secondary_Currents(XYZ,Js,R,ax[0])
+
+        elif Figure1b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[0])
+
+
+    if Figure1a== 'Configuration':
+        Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+        ax[1] = Plot_Primary_Potential(XYZ,Vp,R,ax[1])
+        print 'While figure1 is plotting Configuration, figure2 plots the primary field'
+
+    elif Figure2a == 'Total':
+    
+        if Figure2b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Total_Potential(XYZ,Vt,R,ax[1])
+
+        elif Figure2b == 'ElectricField':
+            ax[0] = Plot_Total_ElectricField(XYZ,Et,R,ax[1])
+
+        elif Figure2b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Total_Currents(XYZ,Jt,R,ax[1])
+
+        elif Figure2b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[1])
+
+
+    elif Figure2a == 'Secondary':
+        
+        if Figure2b == 'Potential':
+            Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+            ax[0] = Plot_Secondary_Potential(XYZ,Vs,R,ax[1])
+
+        elif Figure2b == 'ElectricField':
+            ax[0] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[1])
+
+        elif Figure2b == 'CurrentDensity':
+            Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+            ax[0] = Plot_Secondary_Currents(XYZ,Js,R,ax[1])
+
+        elif Figure2b == 'ChargesDensity':
+            rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+            ax[0] = Plot_ChargesDensity(XYZ,rho,R,ax[1])
+
+    plt.tight_layout(True)
+    plt.show()
+
+#Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey
+def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example):
+
+    sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
+    sig2 = 10.**log_sig2
+    E0   = 1.           # inducing field strength in V/m
+    n = 100             #level of discretisation
+    xr = np.linspace(-200., 200., n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+    PlotOpt = 'Total'
+
+    if matching_spheres_example:
+        sig0 = 10.**(-3)
+        sig1 = 10.**(-2)
+        sig2 = 1.310344828 * 10**(-3)
+        R0   = 20.
+        R1   = 40.
+
+        two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
+
+    else:
+        two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
+
+    plt.tight_layout(True)
+    plt.show()
+
+def run(plotIt=True):
+    sig0 = -3.          # conductivity of the wholespace
+    sig1 = -1.         # conductivity of the sphere
+    sig0, sig1 = conductivity_log_wrapper(sig0,sig1)
+    R    = 50.          # radius of the sphere
+    E0   = 1.           # inducing field strength
+    n = 100             #level of discretisation
+    xr = np.linspace(-2.*R, 2.*R, n) # X-axis discretization
+    yr = xr.copy()      # Y-axis discretization
+    zr = np.r_[0]          # identical to saying `zr = np.array([0])`
+    XYZ = ndgrid(xr,yr,zr) # Space Definition
+            
+    Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
+    Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
+    Jt,Jp,Js, = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
+    rho = get_ChargesDensity(XYZ,sig0,sig1,R,Ep)
+
+    if plotIt:
+        fig, ax = plt.subplots(2,5,figsize=(50,10))
+        ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1])
+        ax[1,0] = Plot_Primary_Potential(XYZ,Vp,R,ax[1,0])
+        ax[0,1] = Plot_Total_Potential(XYZ,Vt,R,ax[0,1])
+        ax[1,1] = Plot_Secondary_Potential(XYZ,Vs,R,ax[1,1])
+        ax[0,2] = Plot_Total_ElectricField(XYZ,Et,R,ax[0,2])
+        ax[1,2] = Plot_Secondary_ElectricField(XYZ,Es,R,ax[1,2])
+        ax[0,3] = Plot_Total_Currents(XYZ,Jt,R,ax[0,3])
+        ax[1,3] = Plot_Secondary_Currents(XYZ,Js,R,ax[1,3])
+        ax[0,4] = Plot_Primary_Potential(XYZ,Vp,R,ax[0,4])
+        ax[1,4] = Plot_ChargesDensity(XYZ,rho,R,ax[1,4])
+
+
+        plt.show()
+
+if __name__ == '__main__':
+    run()
+
+

SimPEG/MT/Utils/ediFilesUtils.py

@@ -7,17 +7,16 @@ from SimPEG.MT.Utils.dataUtils import rec2ndarr
 # Import modules
 import numpy as np
 import os, sys, re
-try:
-    import osr
-except ImportError as e:
-    print 'Could not import osr, missing the gdal package'
-    pass
+
 
 class EDIimporter:
     """
     A class to import EDIfiles.
 
     """
+
+
+    # Define data converters
     _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
     _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
 
@@ -26,8 +25,8 @@ class EDIimporter:
     comps = None
 
     # Hidden properties
-    _outEPSG = None
-    _2out = None
+    _outEPSG = None # Project info
+    _2out = None    # The projection operator
 
 
     def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -113,6 +112,12 @@ class EDIimporter:
     # nOutData=length(obj.data);
     # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
     def _transfromPoints(self,longD,latD):
+        # Import the coordinate projections
+        try:
+            import osr
+        except ImportError as e:
+            print 'Could not import osr, missing the gdal package\nCan not project coordinates'
+            raise e
         # Coordinates convertor
         if self._2out is None:
             src = osr.SpatialReference()

SimPEG/Mesh/CylMesh.py

@@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                 raise NotImplementedError('wrapping in the averaging is not yet implemented')
         return self._aveF2CCV
 
-    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
         """
             Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
         """
@@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
 
 
+        if locTypeTo is None:
+            locTypeTo = locType
+
         if locType == 'F':
             # do this three times for each component
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
-            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
             return sp.vstack((X,Y,Z))
         if locType == 'E':
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
-            Z = spzeros(Mrect.nEz, self.nE)
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
+            Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
             return sp.vstack((X,Y,Z))
 
-        grid = getattr(Mrect, 'grid' + locType)
+        grid = getattr(Mrect, 'grid' + locTypeTo)
         # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
         theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
         theta[theta < 0] += np.pi*2.0
@@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                         'Ex': Mrect.tangents[:Mrect.nEx,:],
                         'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
                         'Ez': Mrect.tangents[-Mrect.nEz:,:],
-                    }[locType]
+                    }[locTypeTo]
             if 'F' in locType:
                 normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
                 proj = ( normals * dotMe ).sum(axis=1)

SimPEG/Mesh/MeshIO.py

@@ -21,10 +21,9 @@ class TensorMeshIO(object):
                 if '*' in seg:
                     st = seg
                     sp = seg.split('*')
-                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
+                    re = int(sp[0])*(' ' + sp[1])
                     line = line.replace(st,re.strip())
             return np.array(line.split(),dtype=float)
-
         # Read the file as line strings, remove lines with comment = !
         msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
 

SimPEG/Optimization.py

@@ -1003,7 +1003,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 # perturb inactive set off of bounds so that they are included in the step
                 delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
 
-
             # 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.

SimPEG/Regularization.py

@@ -311,6 +311,9 @@ class BaseRegularization(object):
             tmp = indActive
             indActive = np.zeros(mesh.nC, dtype=bool)
             indActive[tmp] = True
+        if indActive is not None and mapping is None:
+            mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
         self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
@@ -728,14 +731,14 @@ class Sparse(Simple):
     @property
     def W(self):
         """Full regularization matrix W"""
-        #if getattr(self, '_W', None) is None:
-        wlist = (self.Wsmall, self.Wsmooth)
-        #self._W = sp.vstack(wlist)
-        return sp.vstack(wlist)
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Wsmall, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
 
     def R(self, f_m , eps, exponent):
 
-        eta = (eps**(1-exponent/2.))**0.5
-        r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
+        eta = (eps**(1.-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
 
         return r

SimPEG/Utils/__init__.py

@@ -7,3 +7,4 @@ from CounterUtils import *
 import ModelBuilder
 import SolverUtils
 from coordutils import *
+from plottingUtils import *

SimPEG/Utils/io_utils.py

@@ -0,0 +1,137 @@
+from SimPEG import np, Mesh
+import time as tm
+import vtk, vtk.util.numpy_support as npsup
+import re
+
+def read_GOCAD_ts(tsfile):
+    """
+
+    Read GOCAD triangulated surface (*.ts) file
+    INPUT:
+    tsfile: Triangulated surface
+
+    OUTPUT:
+    vrts : Array of vertices in XYZ coordinates [n x 3]
+    trgl : Array of index for triangles [m x 3]. The order of the vertices
+            is important and describes the normal
+            n = cross( (P2 - P1 ) , (P3 - P1) )
+
+    Author: @fourndo
+
+
+    .. note::
+
+        Remove all attributes from the GoCAD surface before exporting it!
+
+    """
+
+
+    fid = open(tsfile,'r')
+    line = fid.readline()
+
+    # Skip all the lines until the vertices
+    while re.match('TFACE',line)==None:
+        line = fid.readline()
+
+    line = fid.readline()
+    vrtx = []
+
+    # Run down all the vertices and save in array
+    while re.match('VRTX',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[2:5])
+        vrtx.append(temp.astype(np.float))
+
+        # Read next line
+        line = fid.readline()
+
+    vrtx = np.asarray(vrtx)
+
+    # Skip lines to the triangles
+    while re.match('TRGL',line)==None:
+        line = fid.readline()
+
+    # Run down the list of triangles
+    trgl = []
+
+    # Run down all the vertices and save in array
+    while re.match('TRGL',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[1:4])
+        trgl.append(temp.astype(np.int))
+
+        # Read next line
+        line = fid.readline()
+
+    trgl = np.asarray(trgl)
+
+    return vrtx, trgl
+
+def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True):
+    """"
+    Function to read gocad polystructure file and output indexes of mesh with in the structure.
+
+    """
+    # Adjust the index
+    trgl = trgl - 1
+
+    # Make vtk pts
+    ptsvtk = vtk.vtkPoints()
+    ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1))
+
+    # Make the polygon connection
+    polys = vtk.vtkCellArray()
+    for face in trgl:
+        poly = vtk.vtkPolygon()
+        poly.GetPointIds().SetNumberOfIds(len(face))
+        for nrv, vert in enumerate(face):
+            poly.GetPointIds().SetId(nrv,vert)
+        polys.InsertNextCell(poly)
+
+    # Make the polydata, structure of connections and vrtx
+    polyData = vtk.vtkPolyData()
+    polyData.SetPoints(ptsvtk)
+    polyData.SetPolys(polys)
+
+    # Make implicit func
+    ImpDistFunc = vtk.vtkImplicitPolyDataDistance()
+    ImpDistFunc.SetInput(polyData)
+
+    # Convert the mesh
+    vtkMesh = vtk.vtkRectilinearGrid()
+    vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz)
+    vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1))
+    vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1))
+    vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1))
+    # Add indexes
+    vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1)
+    vtkInd.SetName('Index')
+    vtkMesh.GetCellData().AddArray(vtkInd)
+
+    extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor
+    extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) #
+    extractImpDistRectGridFilt.SetInputData(vtkMesh)
+
+    if boundaries is True:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOff()
+
+    if internal is True:
+        extractImpDistRectGridFilt.ExtractInsideOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractInsideOff()
+
+    print "Extracting indices from grid..."
+    # Executing the pipe
+    extractImpDistRectGridFilt.Update()
+
+    # Get index inside
+    insideGrid = extractImpDistRectGridFilt.GetOutput()
+    insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))
+
+
+    # Return the indexes inside
+    return insideGrid

SimPEG/Utils/plottingUtils.py

@@ -0,0 +1,3 @@
+# Plot Tree!
+# Plot SphereSetup
+# Plot LayerEarth

docs/examples/DC_Forward_PseudoSection.rst

@@ -22,7 +22,7 @@ radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
 stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
 dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
-Created by @fourndo on Mon Feb 01 19:28:06 2016
+Created by @fourndo
 
 
 

docs/examples/MT_1D_analytic_nlayer_Earth.rst

@@ -0,0 +1,21 @@
+.. _examples_MT_1D_analytic_nlayer_Earth:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+MT 1D analytic nlayer Earth
+===========================
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.MT_1D_analytic_nlayer_Earth.run()
+
+.. literalinclude:: ../../SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+    :language: python
+    :linenos:

docs/examples/sphereElectrostatic_example.rst

@@ -0,0 +1,21 @@
+.. _examples_sphereElectrostatic_example:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+sphereElectrostatic example
+===========================
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.sphereElectrostatic_example.run()
+
+.. literalinclude:: ../../SimPEG/Examples/sphereElectrostatic_example.py
+    :language: python
+    :linenos:

tests/base/test_regularization.py

@@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase):
                     elif mesh.dim == 3:
                         indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
-
                     for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
-                        reg = r(mesh, mapping=mapping, indActive=indAct)
+                        reg = r(mesh, indActive=indAct)
                         m = np.random.rand(mesh.nC)[indAct]
                         reg.mref = np.ones_like(m)*np.mean(m)
 

tests/mesh/test_cylMesh.py

@@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase):
 
         assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
 
+    def test_getInterpMatCartMesh_Cells2Nodes(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
     def test_getInterpMatCartMesh_Faces(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase):
         assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
+    def test_getInterpMatCartMesh_Faces2Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E')
+        mf = np.ones(Mc.nF)
+
+        ecart = Pf2e * mf
+
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 1) < TOL
+        assert np.abs(float(excc[indY]) - 0) < TOL
+        assert np.abs(float(eycc[indX]) - 0) < TOL
+        assert np.abs(float(eycc[indY]) - 1) < TOL
+        assert np.abs((ezcc - 1).sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
     def test_getInterpMatCartMesh_Edges(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase):
         Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
         me = np.ones(Mc.nE)
 
-        erect = Pe * me
+        ecart = Pe * me
 
-        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
-        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
-        ezcc = Mr.r(erect, 'E', 'Ez')
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges2Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F')
+        me = np.ones(Mc.nE)
+
+        frect = Pe2f * me
+
+        excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        ezcc = Mr.r(frect, 'F', 'Fz')
 
         indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
         indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])

tests/mt/test_Problem3D_againstAnalytic.py

@@ -242,9 +242,6 @@ class TestAnalytics(unittest.TestCase):
     def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
     def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
 
-    # Do a derivative test
-    def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
-
     # Do a derivative test of Jvec
     # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
     # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))