From f3c96261333287b1853ff94ac9e6f7725cafe19a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 6 Apr 2016 17:35:19 -0700 Subject: [PATCH] placeholder for plotting utils, leverage a bit more simpeg functionality in examples --- .../Examples/MT_1D_analytic_nlayer_Earth.py | 181 +++++----- .../Examples/sphereElectrostatic_example.py | 323 +++++++++--------- SimPEG/Utils/__init__.py | 1 + SimPEG/Utils/plottingUtils.py | 3 + 4 files changed, 251 insertions(+), 257 deletions(-) create mode 100644 SimPEG/Utils/plottingUtils.py diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py index 0cccc5f1..62ab2a68 100644 --- a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py +++ b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py @@ -2,7 +2,7 @@ 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 +from SimPEG.EM.Utils import k, omega """ MT1D: n layered earth problem @@ -16,51 +16,45 @@ This code compute the analytic response of a n-layered Earth to a plane wave (Ma 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 +\\\(\\\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 +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: +The iteration from one layer to another is ensure by: -\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = +\\(\\ \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 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]), +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.]), +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.]), +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.]), +eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]), np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1) ,decimals=1)) @@ -70,17 +64,8 @@ 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 +top = lambda thick: np.cumsum(thick) #Propagation Matrix and theirs inverses @@ -105,36 +90,36 @@ 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' , '.' , '*' ] - + 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]: @@ -144,39 +129,39 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): 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.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.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', @@ -195,7 +180,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): 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', @@ -232,7 +217,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) ) - + ax.invert_yaxis() return ax @@ -240,83 +225,83 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): #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) + + 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): + + 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) +#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): @@ -335,44 +320,44 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): 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)))), + 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)))), + 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()]) @@ -380,12 +365,12 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): 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.]) @@ -395,7 +380,7 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2 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 @@ -410,9 +395,9 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,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 @@ -430,14 +415,14 @@ def run(n,plotIt=True): if plotIt: - PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True) + 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/sphereElectrostatic_example.py b/SimPEG/Examples/sphereElectrostatic_example.py index 7ff1ead1..7bd633c5 100644 --- a/SimPEG/Examples/sphereElectrostatic_example.py +++ b/SimPEG/Examples/sphereElectrostatic_example.py @@ -32,7 +32,7 @@ sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0) def conductivity_log_wrapper(log_sig0,log_sig1): sig0 = 10.**log_sig0 sig1 = 10.**log_sig1 - + return sig0,sig1 # Examples @@ -54,7 +54,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): 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 ) @@ -85,7 +85,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): 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.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', @@ -96,7 +96,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): 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.)), @@ -114,7 +114,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): else: ax.set_xticklabels([]) ax.set_yticklabels([]) - ax.text(-0.05,-10,'$\sigma_0$',fontsize=14) + 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) @@ -128,8 +128,8 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): ax.set_aspect('equal') - - + + return ax def get_Conductivity(XYZ,sig0,sig1,R): @@ -138,61 +138,61 @@ def get_Conductivity(XYZ,sig0,sig1,R): ''' 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): +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) @@ -205,19 +205,19 @@ def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax): 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) @@ -230,16 +230,16 @@ def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax): 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)) @@ -254,30 +254,30 @@ def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax): 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, + 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); @@ -285,16 +285,16 @@ def get_ElectricField(XYZ,sig0,sig1,R,E0): 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] @@ -302,7 +302,7 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): 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) @@ -310,22 +310,22 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): 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 + +#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] @@ -333,7 +333,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): 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) @@ -341,7 +341,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): 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') @@ -349,53 +349,53 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): 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, + 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[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') @@ -403,30 +403,30 @@ def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax): 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') @@ -434,52 +434,52 @@ def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax): 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, + 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') @@ -492,11 +492,11 @@ def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax): 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 @@ -513,20 +513,20 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): #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 = 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)) + 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) @@ -535,36 +535,36 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): 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 + + #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) @@ -580,10 +580,10 @@ def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend #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) @@ -594,8 +594,8 @@ def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend 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) + + #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' ) @@ -608,36 +608,36 @@ def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend #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') + 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') + + 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 @@ -648,45 +648,45 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig fig, ax = plt.subplots(1,2,figsize=(18,6)) - #Setup figure 1 with options Configuration, Total or Secondary, + #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': + + elif Figure2a == 'Total': if Figure2b == 'Potential': ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1]) @@ -699,8 +699,8 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig elif Figure2b == 'ChargesDensity': ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2a == 'Secondary': + + elif Figure2a == 'Secondary': if Figure2b == 'Potential': ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) @@ -715,10 +715,10 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig 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 @@ -728,12 +728,12 @@ def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,x 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) + sig0 = 10.**(-3) + sig1 = 10.**(-2) sig2 = 1.310344828 * 10**(-3) - R0 = 20. + R0 = 20. R1 = 40. two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt) @@ -744,9 +744,7 @@ def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,x plt.tight_layout(True) plt.show() - - -if __name__ == '__main__': +def run(plotIt=True): sig0 = -3. # conductivity of the wholespace sig1 = -1. # conductivity of the sphere sig0, sig1 = conductivity_log_wrapper(sig0,sig1) @@ -758,18 +756,25 @@ if __name__ == '__main__': 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]) - + 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,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]) + else: + get_Potential(XYZ,sig0,sig1,R,E0) # This is so travis tests it + + + + plt.show() + +if __name__ == '__main__': - plt.show() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 18c1994f..3f5d62c8 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,3 +7,4 @@ from CounterUtils import * import ModelBuilder import SolverUtils from coordutils import * +from plottingUtils import * diff --git a/SimPEG/Utils/plottingUtils.py b/SimPEG/Utils/plottingUtils.py new file mode 100644 index 00000000..72ab7654 --- /dev/null +++ b/SimPEG/Utils/plottingUtils.py @@ -0,0 +1,3 @@ +# Plot Tree! +# Plot SphereSetup +# Plot LayerEarth