wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-28 01:00:33 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
7f77cc2ea350d7d3f75941b723ab0b78db799679
simpeg/SimPEG/Examples/sphereElectrostatic_example.py
T
Thibaut Astic 7f77cc2ea3 add link to sphere webpage
2016-02-16 10:49:13 -08:00

797 lines
30 KiB
Python
Raw Blame History

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.
Please visit http://em.geosci.xyz/en/latest/content/maxwell2_steady_state/electrostatic_sphere.html
for more examples using this code.
'''
# 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,ax):
#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
ax[0] = get_Setup(XYZ,sig0,sig1,R0,E0,ax[0],True,[0.6,0.1,0.1])
#Plotting the Configuration 1
ax[1] = get_Setup(XYZ,sig0,sig2,R1,E0,ax[1],True,[0.1,0.1,0.6])
#Plotting the Data (Legends)
ax[2].set_title('Potential Differences',fontsize=ftsize_title)
ax[2].set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
ax[2].set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
ax[2].tick_params(labelsize=ftsize_axis)
ax[2].grid()
if PlotOpt == 'Total':
ax[3]= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax[3])
ax[4]= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax[4])
#Plot the Data (from Configuration 0)
gphy0 = ax[2].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 = ax[2].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' )
ax[2].legend(('Left Model Response','Right Model Response'),loc=4)
elif PlotOpt == 'Secondary':
#plot the secondary potentials
ax[3]= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax[3])
ax[4]= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax[4])
#Plot the data(from configuration 0)
gphy0 = ax[2].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 = ax[2].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' )
ax[2].legend(('Left Model Response','Right Model Response'),loc=4 )
else:
print('What dont you get? Total or Secondary?')
#Legends
ax[3].plot(MP0[:,0],MP0[:,1],color='gray')
Dip_Midpoint0 = ax[3].scatter(MP0[:,0],MP0[:,1],color='black')
Electrodes0 = ax[3].scatter(EL0[:,0],EL0[:,1],color='red')
ax[3].legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
ax[4].plot(MP1[:,0],MP1[:,1],color='gray')
Dip_Midpoint1 = ax[4].scatter(MP1[:,0],MP1[:,1],color='black')
Electrodes1 = ax[4].scatter(EL1[:,0],EL1[:,1],color='red')
ax[4].legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
return ax
#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'
#Initializing the figure
fig = plt.figure(figsize=(20,20))
ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) #Configuration Conductive Sphere
ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) #Configuration Resistive Sphere
ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) # Data
ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) #Potential Conductive Sphere
ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) #Potential Resistive Potential
ax = [ax0,ax1,ax2,ax3,ax4]
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,ax)
else:
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt,ax)
plt.tight_layout(True)
plt.show()
def run(PlotIt=True):
sig0 = -3. # conductivity of the wholespace
sig1 = -1. # conductivity of the sphere
sig0, sig1 = conductivity_log_wrapper(sig0,sig1)
R = 50. # radius of the sphere
E0 = 1. # inducing field strength
n = 100 #level of discretisation
xr = np.linspace(-2.*R, 2.*R, n) # X-axis discretization
yr = xr.copy() # Y-axis discretization
zr = np.r_[0] # identical to saying `zr = np.array([0])`
XYZ = ndgrid(xr,yr,zr) # Space Definition
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep)
if PlotIt:
fig, ax = plt.subplots(2,5,figsize=(50,10))
ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1])
ax[1,0] = Plot_Primary_Potential(XYZ,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()
return Vt,Vp,Vs,Et,Ep,Es,Jt,Jp,Js,rho
if __name__ == '__main__':
run()
Reference in New Issue View Git Blame Copy Permalink