diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index acccb224..8f93c70a 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -242,10 +242,10 @@ 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) ) @@ -266,11 +266,11 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # 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): @@ -292,12 +292,12 @@ class update_IRLS(InversionDirective): gamma = None phi_m_last = None phi_d_last = None - + def initialize(self): # Scale the regularization for changes in norm if getattr(self, 'phi_m_last', None) is not None: - + self.reg.curModel = self.invProb.curModel self.reg.gamma = 1. phim_new = self.reg.eval(self.invProb.curModel) @@ -305,10 +305,10 @@ class update_IRLS(InversionDirective): self.reg.curModel = self.invProb.curModel self.reg.gamma = self.gamma - + if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d - + def endIter(self): # Cool the threshold parameter if getattr(self, 'factor', None) is not None: @@ -321,35 +321,35 @@ class update_IRLS(InversionDirective): # Get phi_m at the end of current iteration self.phi_m_last = self.invProb.phi_m_last - + # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - + # Temporarely set gamma to 1. self.reg.gamma = 1. - + # Compute change in model objective function and update scaling phim_new = self.reg.eval(self.invProb.curModel) self.reg.gamma = self.phi_m_last / phim_new - - - # TO DO: Check optimization class, data misfit not matching reality + + + # TO DO: Check optimization class, data misfit not matching reality #dpred = self.prob.fields(self.invProb.curModel) #phid = self.invProb.dmisfit.eval(self.invProb.curModel) #print self.survey.std[0] #print phid #print self.invProb.phi_d #print self.invProb.phi_d_last - + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d class update_lin_PreCond(InversionDirective): - - + + def endIter(self): # Cool the threshold parameter - + 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. @@ -362,20 +362,20 @@ class update_Wj(InversionDirective): Create approx-sensitivity base weighting """ k = None - + def endIter(self): - + if self.opt.iter == 2: m = self.invProb.curModel if self.k is None: self.k = int(self.survey.nD/10) - + def JtJv(v): Jv = self.prob.Jvec(m, v) - + return self.prob.Jtvec(m,Jv) - + JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k) JtJdiag = JtJdiag / max(JtJdiag) diff --git a/SimPEG/Examples/EM_FDEM_SusEffects.py b/SimPEG/Examples/EM_FDEM_SusEffects.py deleted file mode 100644 index 1abbb16f..00000000 --- a/SimPEG/Examples/EM_FDEM_SusEffects.py +++ /dev/null @@ -1,148 +0,0 @@ -from SimPEG import * -from SimPEG import EM -from pymatsolver import MumpsSolver -from scipy.constants import mu_0 - -def run(plotIt=True): - """ - EM: FDEM: Effects of susceptibility - =================================== - - When airborne freqeuncy domain EM (AFEM) survey is flown over - the earth including significantly susceptible bodies (magnetite-rich rocks), - negative data is often observed in the real part of the lowest frequency - (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization - occurs due to a susceptible body when the magnetic field is applied. - - To clarify what is happening in the earth when we are exciting the earth with - a loop source in the frequency domain we run three forward modelling: - - - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility - - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity - - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field - - We plot vector magnetic fields in the earth. For secondary fields we provide - F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following - figure show both real and parts. - - """ - # Generate Cylindrical mesh - cs, ncx, ncz, npad = 5, 25, 24, 20. - hx = [(cs,ncx), (cs,npad,1.3)] - hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] - mesh = Mesh.CylMesh([hx,1,hz], '00C') - sighalf = 1e-3 - sigma = np.ones(mesh.nC)*1e-8 - sigmahomo = sigma.copy() - mu = np.ones(mesh.nC)*mu_0 - sigma[mesh.gridCC[:,-1]<0.] = sighalf - blkind = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50)) - sigma[blkind] = 1e-1 - mu[blkind] = mu_0*1.1 - offset = 0. - frequency = np.r_[10., 100., 1000.] - rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr') - rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi') - srcLists = [] - nfreq = frequency.size - for ifreq in range(nfreq): - src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.) - srcLists.append(src) - survey = EM.FDEM.Survey(srcLists) - iMap = Maps.IdentityMap(nP=int(mesh.nC)) - # Use PhysPropMap - maps = [('sigma', iMap), ('mu', iMap)] - prob = EM.FDEM.Problem_b(mesh, mapping=maps) - prob.Solver = MumpsSolver - survey.pair(prob) - m = np.r_[sigma, mu] - survey0 = EM.FDEM.Survey(srcLists) - prob0 = EM.FDEM.Problem_b(mesh, mapping=maps) - prob0.Solver = MumpsSolver - survey0.pair(prob0) - m = np.r_[sigma, mu] - m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0] - m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0] - # Anomalous conductivity and susceptibility - F = prob.fields(m) - # Only anomalous conductivity - F0 = prob.fields(m0) - # Primary field - F00 = prob.fields(m00) - - if plotIt: - import matplotlib.pyplot as plt - def vizfields(ifreq=0, primsec="secondary",realimag="real"): - - titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"] - actind = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200)) - - if primsec=="secondary": - bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') - bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim - bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim - elif primsec=="primary": - bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') - bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') - - XYZ = mesh.gridCC[actind,:] - X = XYZ[:,0].reshape((31,43), order='F') - Z = XYZ[:,2].reshape((31,43), order='F') - bx = bCC[actind,0].reshape((31,43), order='F') - bz = bCC[actind,1].reshape((31,43), order='F') - bx0 = bCC0[actind,0].reshape((31,43), order='F') - bz0 = bCC0[actind,1].reshape((31,43), order='F') - - bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F') - bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F') - - absbreal = np.sqrt(bx.real**2+bz.real**2) - absbimag = np.sqrt(bx.imag**2+bz.imag**2) - absb0real = np.sqrt(bx0.real**2+bz0.real**2) - absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2) - - absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2) - absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2) - - fig = plt.figure(figsize=(15,5)) - ax1 = plt.subplot(131) - ax2 = plt.subplot(132) - ax3 = plt.subplot(133) - typefield="real" - scale=20 - if realimag=="real": - ax1.contourf(X, Z,np.log10(absbreal), 100) - ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5) - ax2.contourf(X, Z,np.log10(absb0real), 100) - ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5) - ax3.contourf(X, Z,np.log10(absbrealsec), 100) - ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5) - elif realimag=="imag": - ax1.contourf(X, Z,np.log10(absbimag), 100) - ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5) - ax2.contourf(X, Z,np.log10(absb0imag), 100) - ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5) - ax3.contourf(X, Z,np.log10(absbimagsec), 100) - ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5) - - ax = [ax1, ax2, ax3] - ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18) - ax2.text(30, 50, primsec, color="k", fontsize=18) - ax1.text(30, 50, realimag, color="k", fontsize=18) - for i, axtemp in enumerate(ax): - axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3) - - axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3) - axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3) - axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3) - axtemp.set_ylim(-200, 100.) - axtemp.set_xlim(10, 100.) - axtemp.set_title(titles[i]) - plt.show() - return fig, ax - fig1, ax1 = vizfields(1, primsec="primary", realimag="real") - fig2, ax2 = vizfields(1, primsec="secondary", realimag="real") - fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag") - -if __name__ == '__main__': - run() diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py deleted file mode 100644 index 0cccc5f1..00000000 --- a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py +++ /dev/null @@ -1,443 +0,0 @@ -from scipy.constants import epsilon_0, mu_0 -import matplotlib.pyplot as plt -import numpy as np -from ipywidgets import * -#from SimPEG.EM.Utils import k, omega - -""" -MT1D: n layered earth problem -***************************** - -Author: Thibaut Astic -Contact: thast@eos.ubc.ca -Date: January 2016 - -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. - -""" - -#Frequency conversion -omega = lambda f: 2.*np.pi*f - -#Evaluate k wavenumber -k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f)) - -#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 Perties 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 -def top(thick): - topv= np.zeros(len(thick)+1) - - topv[0]=-thick[0] - - for i in range(1,len(topv),1): - topv[i] = topv[i-1] + thick[i-1] - - return topv - -#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(mu,sigcm,eps,f) - 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(n,plotIt=True): - # 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(3) - - - - - diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index a95bc194..b06bb3a7 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -6,7 +6,6 @@ import DC_Forward_PseudoSection import DC_PseudoSection_Simulation import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace -import EM_FDEM_SusEffects import EM_Schenkel_Morrison_Casing import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 @@ -20,12 +19,10 @@ 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", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "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"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_PseudoSection_Simulation", "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"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/Examples/sphereElectrostatic_example.py b/SimPEG/Examples/sphereElectrostatic_example.py deleted file mode 100644 index 7ff1ead1..00000000 --- a/SimPEG/Examples/sphereElectrostatic_example.py +++ /dev/null @@ -1,775 +0,0 @@ -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 -Date: December 2015 - -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,sig0,sig1,R,E0,ax): - - Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - - 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,sig0,sig1,R,E0,ax): - - Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - - 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,sig0,sig1,R,E0,ax): - - Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - - 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,sig0,sig1,R,E0,ax): - - Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - - 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,sig0,sig1,R,E0,ax): - - Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - - 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,sig0,sig1,R,E0,ax): - - Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) - - 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,sig0,sig1,R,E0,ax): - - Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) - - 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+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,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,Et,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,sig0,sig1,R,E0,ax): - - xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) - xcirc = xr[np.abs(xr) <= R] - - Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep) - - 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() - - if PlotOpt == 'Total': - ax3= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax3) - ax4= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,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,sig0,sig1,R0,E0,ax3) - ax4= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax4) - - #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 - - 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': - ax[0] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'ElectricField': - ax[0] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'CurrentDensity': - ax[0] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'ChargesDensity': - ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1a == 'Secondary': - - if Figure1b == 'Potential': - ax[0] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'ElectricField': - ax[0] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'CurrentDensity': - ax[0] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[0]) - - elif Figure1b == 'ChargesDensity': - ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) - - - if Figure1a== 'Configuration': - ax[1] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) - print 'While figure1 is plotting Configuration, figure2 plots the primary field' - - elif Figure2a == 'Total': - if Figure2b == 'Potential': - ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'ElectricField': - ax[1] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'CurrentDensity': - ax[1]=Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'ChargesDensity': - ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1]) - - - elif Figure2a == 'Secondary': - if Figure2b == 'Potential': - ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'ElectricField': - ax[1] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'CurrentDensity': - ax[1] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2b == 'ChargesDensity': - ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,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() - - - -if __name__ == '__main__': - 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 - - 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,sig0,sig1,R,E0,ax[1,0]) - ax[0,1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0,1]) - ax[1,1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1,1]) - ax[0,2] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0,2]) - ax[1,2] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1,2]) - ax[0,3] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0,3]) - ax[1,3] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1,3]) - ax[0,4] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[0,4]) - ax[1,4] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1,4]) - - - plt.show() - diff --git a/docs/examples/EM_FDEM_SusEffects.rst b/docs/examples/EM_FDEM_SusEffects.rst deleted file mode 100644 index 3e7dea80..00000000 --- a/docs/examples/EM_FDEM_SusEffects.rst +++ /dev/null @@ -1,41 +0,0 @@ -.. _examples_EM_FDEM_SusEffects: - -.. --------------------------------- .. -.. .. -.. THIS FILE IS AUTO GENEREATED .. -.. .. -.. SimPEG/Examples/__init__.py .. -.. .. -.. --------------------------------- .. - - -EM: FDEM: Effects of susceptibility -=================================== - -When airborne freqeuncy domain EM (AFEM) survey is flown over -the earth including significantly susceptible bodies (magnetite-rich rocks), -negative data is often observed in the real part of the lowest frequency -(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization -occurs due to a susceptible body when the magnetic field is applied. - -To clarify what is happening in the earth when we are exciting the earth with -a loop source in the frequency domain we run three forward modelling: - - - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility - - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity - - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field - -We plot vector magnetic fields in the earth. For secondary fields we provide -F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following -figure show both real and parts. - - - -.. plot:: - - from SimPEG import Examples - Examples.EM_FDEM_SusEffects.run() - -.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py - :language: python - :linenos: diff --git a/docs/examples/MT_1D_analytic_nlayer_Earth.rst b/docs/examples/MT_1D_analytic_nlayer_Earth.rst deleted file mode 100644 index 834121a9..00000000 --- a/docs/examples/MT_1D_analytic_nlayer_Earth.rst +++ /dev/null @@ -1,21 +0,0 @@ -.. _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: