mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-17 11:32:59 +08:00
seogi_macbook
f944f9b76b
2. Add Analytic tests 3. Fix simple bug in PlotSlice for nodal variable 4. Add more analytic function (sphere)
114 lines
3.5 KiB
Python
114 lines
3.5 KiB
Python
import numpy as np
|
|
from scipy.constants import mu_0, pi
|
|
|
|
def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
|
|
"""
|
|
Analytic solution for electric potential from a postive pole
|
|
|
|
Input variables:
|
|
|
|
txloc = a xyz location of A (+) electrode (np.r_[xa, ya, za])
|
|
|
|
rxlocs = [M, N]
|
|
M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
|
|
N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
|
|
|
|
sigma = conductivity (either float or complex)
|
|
flag = "wholsespace" or "halfspace"
|
|
|
|
"""
|
|
M = rxlocs[0]
|
|
N = rxlocs[1]
|
|
|
|
rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
|
|
rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
|
|
|
|
phiM = 1./(4*np.pi*rM*sigma)
|
|
phiN = 1./(4*np.pi*rN*sigma)
|
|
phi = phiM - phiN
|
|
|
|
if flag == "halfspace":
|
|
phi *= 2
|
|
|
|
return phi
|
|
|
|
deg2rad = lambda deg: deg/180.*np.pi
|
|
rad2deg = lambda rad: rad*180./np.pi
|
|
|
|
def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
|
|
flag = "sec", order=12):
|
|
# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
|
|
# flag = "sec", order=12):
|
|
"""
|
|
|
|
Parameters:
|
|
|
|
txloc (array) : current electrode location (x,y,z)
|
|
xc (float) : x center of depressed sphere
|
|
rxloc (array) : electrode locations
|
|
(Nx3 array, # of electrodes)
|
|
radius (float): radius of the sphere (m)
|
|
rho (float) : resistivity of the background (ohm-m)
|
|
rho1 (float) : resistivity of the sphere
|
|
flag (string) : "sec", "total", "prim"
|
|
(default="sec")
|
|
"sec": secondary potential only due to sphere
|
|
"prim": primary potential from the point source
|
|
"total": "sec"+"prim"
|
|
order (float) : maximum order of Legendre polynomial
|
|
(default=12)
|
|
|
|
Written by Seogi Kang (skang@eos.ubc.ca)
|
|
Ph.D. Candidate of University of British Columbia, Canada
|
|
|
|
"""
|
|
|
|
Pleg = []
|
|
# Compute Legendre Polynomial
|
|
for i in range(order):
|
|
Pleg.append(special.legendre(i, monic=0))
|
|
|
|
|
|
rho = 1./sigma
|
|
rho1 = 1./sigma1
|
|
|
|
# Center of the sphere should be aligned in txloc in y-direction
|
|
yc = txloc[1]
|
|
xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]]
|
|
r = np.sqrt( (xyz**2).sum(axis=1) )
|
|
|
|
x0 = abs(txloc[0]-xc)
|
|
|
|
costheta = xyz[:,0]/r * (txloc[0]-xc)/x0
|
|
phi = np.zeros_like(r)
|
|
R = (r**2+x0**2.-2.*r*x0*costheta)**0.5
|
|
# primary potential in a whole space
|
|
prim = rho*1./(4*np.pi*R)
|
|
|
|
if flag =="prim":
|
|
return prim
|
|
|
|
sphind = r < radius
|
|
out = np.zeros_like(r)
|
|
for n in range(order):
|
|
An, Bn = AnBnfun(n, radius, x0, rho, rho1)
|
|
dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind])
|
|
out[~sphind] += dumout
|
|
dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind])
|
|
out[sphind] += dumin
|
|
|
|
out[~sphind] += prim[~sphind]
|
|
|
|
if flag == "sec":
|
|
return out-prim
|
|
elif flag == "total":
|
|
return out
|
|
|
|
def AnBnfun(n, radius, x0, rho, rho1, I=1.):
|
|
const = I*rho/(4*np.pi)
|
|
bunmo = n*rho + (n+1)*rho1
|
|
An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \
|
|
(rho1-rho) / bunmo
|
|
Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo
|
|
return An, Bn
|