mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 04:19:45 +08:00
Fairly major source refactor:
- removed the folder Sources and put those routines inside of SrcUtils. - Broke apart calls for MagDipole, MagDipole_B, CircularLoop
This commit is contained in:
Lindsey Heagy
+171
-98
@@ -1,5 +1,5 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils import SrcUtils
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
|
||||
@@ -91,134 +91,207 @@ class RxFDEM(Survey.BaseRx):
|
||||
# Sources
|
||||
####################################################
|
||||
|
||||
# class SrcFDEM(Survey.BaseSrc):
|
||||
# freq = None
|
||||
# rxPair = RxFDEM
|
||||
# knownSrcTypes = {}
|
||||
|
||||
|
||||
class SrcFDEM(Survey.BaseSrc):
|
||||
#TODO: Break these out into Classes of Sources.
|
||||
|
||||
freq = None #: Frequency (float)
|
||||
|
||||
freq = None
|
||||
rxPair = RxFDEM
|
||||
knownSrcTypes = ['Simple', 'MagDipole'] #TODO: Do we want to just classify them by Magnetic, Electric, Both?
|
||||
|
||||
knownSrcTypes = ['VMD', 'VMD_B', 'CircularLoop', 'Simple']
|
||||
|
||||
radius = None
|
||||
class SrcFDEM_Simple_e(SrcFDEM):
|
||||
"""
|
||||
Simple electric source. It is defined by the user provided vector S_e
|
||||
|
||||
def __init__(self, loc, srcType, freq, rxList):
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, S_e, freq, rxList):
|
||||
self.S_e = S_e
|
||||
self.freq = float(freq)
|
||||
Survey.BaseSrc.__init__(self, loc, srcType, rxList)
|
||||
SrcFDEM.__init__(self, None, 'Simple', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
return None, self.S_e
|
||||
|
||||
src = self
|
||||
freq = src.freq
|
||||
solType = prob._fieldType # Hack, should just ask whether j_m, j_g are defined on edges or faces
|
||||
|
||||
if solType == 'e' or solType == 'b':
|
||||
gridEJx = prob.mesh.gridEx
|
||||
gridEJy = prob.mesh.gridEy
|
||||
gridEJz = prob.mesh.gridEz
|
||||
nEJ = prob.mesh.nE
|
||||
|
||||
gridBHx = prob.mesh.gridFx
|
||||
gridBHy = prob.mesh.gridFy
|
||||
gridBHz = prob.mesh.gridFz
|
||||
nBH = prob.mesh.nF
|
||||
def getSourceDeriv(self, prob, v, adjoint = False):
|
||||
return None, None
|
||||
|
||||
|
||||
class SrcFDEM_Simple_m(SrcFDEM):
|
||||
"""
|
||||
Simple magnetic source. It is defined by the user provided vector S_m
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, S_m, freq, rxList):
|
||||
self.S_m = S_m
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, None, 'Simple', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
return self.S_m, None
|
||||
|
||||
def getSourceDeriv(self, prob, v, adjoint = False):
|
||||
return None, None
|
||||
|
||||
|
||||
class SrcFDEM_Simple(SrcFDEM):
|
||||
"""
|
||||
Simple source. It is defined by the user provided vectors S_m, S_e
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
def __init__(self, S_m, S_e, freq, rxList):
|
||||
self.S_m = S_m
|
||||
self.S_e = S_e
|
||||
SrcFDEM.__init__(self, None, 'Simple', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
return self.S_m, self.S_e
|
||||
|
||||
def getSourceDeriv(self, prob, v, adjoint=None):
|
||||
return None, None
|
||||
|
||||
|
||||
class SrcFDEM_MagDipole(SrcFDEM):
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
def __init__(self, loc, freq, rxList, orientation='Z'):
|
||||
self.freq = float(freq)
|
||||
self.orientation = orientation
|
||||
SrcFDEM.__init__(self, loc, 'MagDipole', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
mui = prob.MfMui
|
||||
|
||||
elif solType == 'h' or solType == 'j':
|
||||
gridEJx = prob.mesh.gridFx
|
||||
gridEJy = prob.mesh.gridFy
|
||||
gridEJz = prob.mesh.gridFz
|
||||
nEJ = prob.mesh.nF
|
||||
|
||||
gridBHx = prob.mesh.gridEx
|
||||
gridBHy = prob.mesh.gridEy
|
||||
gridBHz = prob.mesh.gridEz
|
||||
nBH = prob.mesh.nE
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
mui = prob.MeMuI
|
||||
|
||||
else:
|
||||
NotImplementedError('Only E or F sources')
|
||||
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
|
||||
if src.srcType == 'VMD':
|
||||
SRC = Sources.MagneticDipoleVectorPotential(src.loc, gridEJy, 'y')
|
||||
elif src.srcType == 'CircularLoop':
|
||||
SRC = Sources.MagneticLoopVectorPotential(src.loc, gridEJy, 'y', src.radius)
|
||||
else:
|
||||
raise NotImplementedError('Only VMD and CircularLoop')
|
||||
|
||||
elif prob.mesh._meshType is 'TENSOR':
|
||||
|
||||
if src.srcType == 'VMD':
|
||||
srcfct = Sources.MagneticDipoleVectorPotential
|
||||
SRCx = srcfct(src.loc, gridEJx, 'x')
|
||||
SRCy = srcfct(src.loc, gridEJy, 'y')
|
||||
SRCz = srcfct(src.loc, gridEJz, 'z')
|
||||
|
||||
elif src.srcType == 'VMD_B':
|
||||
srcfct = Sources.MagneticDipoleFields
|
||||
SRCx = srcfct(src.loc, gridBHx, 'x')
|
||||
SRCy = srcfct(src.loc, gridBHy, 'y')
|
||||
SRCz = srcfct(src.loc, gridBHz, 'z')
|
||||
|
||||
elif src.srcType == 'CircularLoop':
|
||||
srcfct = Sources.MagneticLoopVectorPotential
|
||||
SRCx = srcfct(src.loc, gridEJx, 'x', src.radius)
|
||||
SRCy = srcfct(src.loc, gridEJy, 'y', src.radius)
|
||||
SRCz = srcfct(src.loc, gridEJz, 'z', src.radius)
|
||||
else:
|
||||
|
||||
raise NotImplemented('%s srcType is not implemented' % src.srcType)
|
||||
SRC = np.concatenate((SRCx, SRCy, SRCz))
|
||||
a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y')
|
||||
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
srcfct = SrcUtils.MagneticDipoleVectorPotential
|
||||
ax = srcfct(self.loc, gridX, 'x')
|
||||
ay = srcfct(self.loc, gridY, 'y')
|
||||
az = srcfct(self.loc, gridZ, 'z')
|
||||
a = np.concatenate((ax, ay, az))
|
||||
|
||||
# b-forumlation
|
||||
if src.srcType == 'VMD_B':
|
||||
b_0 = SRC
|
||||
S_m = -1j*omega(self.freq)*C*a
|
||||
|
||||
return S_m, None
|
||||
|
||||
|
||||
def getSourceDeriv(self, prob, v, adjoint=None):
|
||||
return None, None
|
||||
|
||||
|
||||
class MagDipole_Bfield(SrcFDEM):
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
def __init__(self, loc, freq, rxList, orientation='Z'):
|
||||
self.freq = float(freq)
|
||||
self.orientation = orientation
|
||||
SrcFDEM.__init__(self, loc, 'MagDipole', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
srcfct = SrcUtils.MagneticDipoleFields
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
bx = srcfct(self.loc, gridX, 'x')
|
||||
bz = srcfct(self.loc, gridZ, 'z')
|
||||
b = np.concatenate((bx,bz))
|
||||
else:
|
||||
a = SRC
|
||||
b_0 = C*a
|
||||
bx = srcfct(self.loc, gridX, 'x')
|
||||
by = srcfct(self.loc, gridY, 'y')
|
||||
bz = srcfct(self.loc, gridZ, 'z')
|
||||
b = np.concatenate((bx,by,bz))
|
||||
|
||||
return -1j*omega(freq)*b_0, None
|
||||
return -1j*omega(self.freq)*b, None
|
||||
|
||||
class SimpleSrcFDEM_e(SrcFDEM):
|
||||
def getSourceDeriv(self, prob, v, adjoint=None):
|
||||
return None, None
|
||||
|
||||
def __init__(self, vec, freq, rxList):
|
||||
self.vec = vec
|
||||
|
||||
class SrcFDEM_CircularLoop(SrcFDEM):
|
||||
|
||||
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
|
||||
def __init__(self, loc, freq, rxList, orientation='Z', radius = 1.):
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, None, 'Simple', freq, rxList)
|
||||
self.orientation = orientation
|
||||
self.radius = radius
|
||||
SrcFDEM.__init__(self, loc, 'MagDipole', rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
return None, self.vec
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y', self.radius)
|
||||
|
||||
else:
|
||||
srcfct = SrcUtils.MagneticDipoleVectorPotential
|
||||
ax = srcfct(self.loc, gridX, 'x', self.radius)
|
||||
ay = srcfct(self.loc, gridY, 'y', self.radius)
|
||||
az = srcfct(self.loc, gridZ, 'z', self.radius)
|
||||
a = np.concatenate((ax, ay, az))
|
||||
|
||||
return -1j*omega(self.freq)*C*a
|
||||
|
||||
|
||||
class SimpleSrcFDEM_m(SrcFDEM):
|
||||
|
||||
def __init__(self, vec, freq, rxList):
|
||||
self.vec = vec
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, None, 'Simple', freq, rxList)
|
||||
|
||||
def getSource(self, prob):
|
||||
return self.vec, None
|
||||
|
||||
####################################################
|
||||
# Survey
|
||||
####################################################
|
||||
|
||||
class SurveyFDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
|
||||
@@ -1,101 +0,0 @@
|
||||
from SimPEG import *
|
||||
from scipy.special import ellipk, ellipe
|
||||
|
||||
def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius):
|
||||
"""
|
||||
Calculate the vector potential of horizontal circular loop
|
||||
at given locations
|
||||
|
||||
:param numpy.ndarray srcLoc: Location of the source(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray I: Input current of the loop
|
||||
:param numpy.ndarray radius: radius of the loop
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius)
|
||||
return np.concatenate(out)
|
||||
|
||||
if isinstance(obsLoc, Mesh.BaseMesh):
|
||||
mesh = obsLoc
|
||||
assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']"
|
||||
return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius)
|
||||
|
||||
srcLoc = np.atleast_2d(srcLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
|
||||
n = obsLoc.shape[0]
|
||||
nSrc = srcLoc.shape[0]
|
||||
|
||||
if component=='z':
|
||||
A = np.zeros((n, nSrc))
|
||||
if nSrc ==1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
else:
|
||||
|
||||
A = np.zeros((n, nSrc))
|
||||
for i in range (nSrc):
|
||||
x = obsLoc[:, 0] - srcLoc[i, 0]
|
||||
y = obsLoc[:, 1] - srcLoc[i, 1]
|
||||
z = obsLoc[:, 2] - srcLoc[i, 2]
|
||||
r = np.sqrt(x**2 + y**2)
|
||||
m = (4 * radius * r) / ((radius + r)**2 + z**2)
|
||||
m[m > 1.] = 1.
|
||||
# m might be slightly larger than 1 due to rounding errors
|
||||
# but ellipke requires 0 <= m <= 1
|
||||
K = ellipk(m)
|
||||
E = ellipe(m)
|
||||
ind = (r > 0) & (m < 1)
|
||||
# % 1/r singular at r = 0 and K(m) singular at m = 1
|
||||
Aphi = np.zeros(n)
|
||||
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
if component == 'x':
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
|
||||
elif component == 'y':
|
||||
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
if nSrc == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import Mesh
|
||||
import matplotlib.pyplot as plt
|
||||
cs = 20
|
||||
ncx, ncy, ncz = 41, 41, 40
|
||||
hx = np.ones(ncx)*cs
|
||||
hy = np.ones(ncy)*cs
|
||||
hz = np.ones(ncz)*cs
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
|
||||
srcLoc = np.r_[0., 0., 0.]
|
||||
Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200)
|
||||
Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200)
|
||||
Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200)
|
||||
A = np.r_[Ax, Ay, Az]
|
||||
B0 = mesh.edgeCurl*A
|
||||
J0 = mesh.edgeCurl.T*B0
|
||||
|
||||
# mesh.plotImage(A, vType = 'Ex')
|
||||
# mesh.plotImage(A, vType = 'Ey')
|
||||
|
||||
mesh.plotImage(B0, vType = 'Fx')
|
||||
mesh.plotImage(B0, vType = 'Fy')
|
||||
mesh.plotImage(B0, vType = 'Fz')
|
||||
|
||||
# # mesh.plotImage(J0, vType = 'Ex')
|
||||
# mesh.plotImage(J0, vType = 'Ey')
|
||||
# mesh.plotImage(J0, vType = 'Ez')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
@@ -1,3 +0,0 @@
|
||||
from magneticDipole import MagneticDipoleVectorPotential
|
||||
from CircularLoop import MagneticLoopVectorPotential
|
||||
from magneticDipole import MagneticDipoleFields
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Solver, Problem
|
||||
from SimPEG.Problem import BaseTimeProblem
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils import SrcUtils
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
from SimPEG import Utils, Mesh
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils import SrcUtils
|
||||
from BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
|
||||
@@ -35,7 +35,7 @@ def getProblem(fdemType, comp):
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
|
||||
Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole(np.r_[0.,0.,0.], freq, [Rx0])
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
import numpy as np
|
||||
from SimPEG import *
|
||||
from scipy.special import ellipk, ellipe
|
||||
from scipy.constants import mu_0, pi
|
||||
from SimPEG import Mesh
|
||||
|
||||
def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
|
||||
"""
|
||||
@@ -53,6 +53,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
|
||||
def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
|
||||
"""
|
||||
Calculate the vector potential of a set of magnetic dipoles
|
||||
@@ -98,3 +99,104 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
|
||||
if nSrc == 1:
|
||||
return B.flatten()
|
||||
return B
|
||||
|
||||
|
||||
|
||||
def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius):
|
||||
"""
|
||||
Calculate the vector potential of horizontal circular loop
|
||||
at given locations
|
||||
|
||||
:param numpy.ndarray srcLoc: Location of the source(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray I: Input current of the loop
|
||||
:param numpy.ndarray radius: radius of the loop
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius)
|
||||
return np.concatenate(out)
|
||||
|
||||
if isinstance(obsLoc, Mesh.BaseMesh):
|
||||
mesh = obsLoc
|
||||
assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']"
|
||||
return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius)
|
||||
|
||||
srcLoc = np.atleast_2d(srcLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
|
||||
n = obsLoc.shape[0]
|
||||
nSrc = srcLoc.shape[0]
|
||||
|
||||
if component=='z':
|
||||
A = np.zeros((n, nSrc))
|
||||
if nSrc ==1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
else:
|
||||
|
||||
A = np.zeros((n, nSrc))
|
||||
for i in range (nSrc):
|
||||
x = obsLoc[:, 0] - srcLoc[i, 0]
|
||||
y = obsLoc[:, 1] - srcLoc[i, 1]
|
||||
z = obsLoc[:, 2] - srcLoc[i, 2]
|
||||
r = np.sqrt(x**2 + y**2)
|
||||
m = (4 * radius * r) / ((radius + r)**2 + z**2)
|
||||
m[m > 1.] = 1.
|
||||
# m might be slightly larger than 1 due to rounding errors
|
||||
# but ellipke requires 0 <= m <= 1
|
||||
K = ellipk(m)
|
||||
E = ellipe(m)
|
||||
ind = (r > 0) & (m < 1)
|
||||
# % 1/r singular at r = 0 and K(m) singular at m = 1
|
||||
Aphi = np.zeros(n)
|
||||
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
if component == 'x':
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
|
||||
elif component == 'y':
|
||||
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
if nSrc == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import Mesh
|
||||
import matplotlib.pyplot as plt
|
||||
cs = 20
|
||||
ncx, ncy, ncz = 41, 41, 40
|
||||
hx = np.ones(ncx)*cs
|
||||
hy = np.ones(ncy)*cs
|
||||
hz = np.ones(ncz)*cs
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
|
||||
srcLoc = np.r_[0., 0., 0.]
|
||||
Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200)
|
||||
Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200)
|
||||
Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200)
|
||||
A = np.r_[Ax, Ay, Az]
|
||||
B0 = mesh.edgeCurl*A
|
||||
J0 = mesh.edgeCurl.T*B0
|
||||
|
||||
# mesh.plotImage(A, vType = 'Ex')
|
||||
# mesh.plotImage(A, vType = 'Ey')
|
||||
|
||||
mesh.plotImage(B0, vType = 'Fx')
|
||||
mesh.plotImage(B0, vType = 'Fy')
|
||||
mesh.plotImage(B0, vType = 'Fz')
|
||||
|
||||
# # mesh.plotImage(J0, vType = 'Ex')
|
||||
# mesh.plotImage(J0, vType = 'Ey')
|
||||
# mesh.plotImage(J0, vType = 'Ez')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
@@ -1,4 +1,5 @@
|
||||
# import Sources
|
||||
# import Ana
|
||||
# import Solver
|
||||
import EMUtils
|
||||
import EMUtils
|
||||
import SrcUtils
|
||||
@@ -2,6 +2,5 @@
|
||||
import TDEM
|
||||
import FDEM
|
||||
import Base
|
||||
import Sources
|
||||
import Analytics
|
||||
import Utils
|
||||
|
||||
Reference in New Issue
Block a user