wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-28 11:26:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

placeholder for plotting utils, leverage a bit more simpeg functionality in examples

Browse Source
This commit is contained in:
Lindsey Heagy
2016-04-06 17:35:19 -07:00
parent ea0500e056
commit f3c9626133
4 changed files with 251 additions and 257 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+83 -98
View File
@@ -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)
SimPEG/Examples/sphereElectrostatic_example.py
+164 -159
View File
@@ -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()
SimPEG/Utils/__init__.py
+1
View File
@@ -7,3 +7,4 @@ from CounterUtils import *
import ModelBuilder
import SolverUtils
from coordutils import *
from plottingUtils import *
SimPEG/Utils/plottingUtils.py
+3
View File
@@ -0,0 +1,3 @@
# Plot Tree!
# Plot SphereSetup
# Plot LayerEarth