In [1]:

    

from SimPEG import *
from scipy.constants import mu_0


def omega(freq):
    """Change frequency to angular frequency, omega"""
    return 2.*np.pi*freq


    

In [2]:

    

M = Mesh.TensorMesh([[(100.,32)],[(100.,34)],[(100.,18)]], x0='CCC')


    

In [3]:

    

print M


    

    




  ---- 3-D TensorMesh ----  
   x0: -1600.00
   y0: -1700.00
   z0: -900.00
  nCx: 32
  nCy: 34
  nCz: 18
   hx: 32*100.00
   hy: 34*100.00
   hz: 18*100.00

In [4]:

    

# Setup the model
conds = [1e-2,1]
sig = Utils.ModelBuilder.defineBlock(M.gridCC,[0,0,-300],[500,500,-100],conds)
sig[M.gridCC[:,2]>0] = 1e-8
sigBG = np.zeros(M.nC) + conds[0]
sigBG[M.gridCC[:,2]>0] = 1e-8
colorbar(M.plotImage(log10(sig)))


    

    
Out[4]:





<matplotlib.colorbar.Colorbar instance at 0x7f55faf27290>






    

In [5]:

    

# Get the mass matrix 
# The model
Msig = M.getEdgeInnerProduct(sig)
MsigBG = M.getEdgeInnerProduct(sigBG)

Mmu = M.getFaceInnerProduct(mu_0, invProp=True)


    

In [6]:

    

freq = 1e1
C = M.edgeCurl
A = C.T*Mmu*C - 1j*omega(freq)*Msig
ABG = C.T*Mmu*C - 1j*omega(freq)*MsigBG


    

In [50]:

    

# Need to solve x and y polarizations of the source.
from simpegMT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_x = np.zeros(M.vnEx,dtype=complex)
ey_x = np.zeros((M.nEy,1),dtype=complex)
ez_x = np.zeros((M.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in arange(M.vnEx[0]):
    for j in arange(M.vnEx[2]):
        ex_x[i,j,:] = e0_1d
eBG_x = np.vstack((simpeg.Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))
rhs_x = ABG.dot(eBG_x)


    

In [12]:

    

# Setup y (north) polarization (_y)
ex_y = np.zeros((M.vnEx,1))
ey_y = np.zeros(M.vnEy)
ez_y = np.zeros((M.vnEz,1))
# Assign the source to ex_x
for i in arange(M.vnEy[0]):
    for j in arange(M.vnEy[2]):
        ey_y[i,j,:] = e0_1d
eBG_y = np.vstack((ex_y,simpeg.Utils.mkvc(M.r(ey_y,'Ey','Ey','V'),2),ez_y))
rhs_y = ABG.dot(eBG_y)


    

In [ ]:

    

# Solve the systems for each polarization
Ainv = simpeg.SolverCG(A)
e_x = Ainv*rhs_x

e_y = Ainv*rhs_x


    

In [28]:

    




    

In [36]:

    




    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-36-4c5017fec386> in <module>()
      5 for i in arange(M.vnEx[0]):
      6     for j in arange(M.vnEx[2]):
----> 7         ex_x[i,j,:] = e0_1d

NameError: name 'e0_1d' is not defined

In [34]:

    

# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh(M.hz,M.x0[2])
e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_x = np.zeros(M.vnEx)
ey_x = np.zeros(M.vnEy)
ez_x = np.zeros(M.vnEz)
# Assign the source to ex_x
for i in arange(M.vnEx[0]):
    for j in arange(M.vnEx[2]):
        ex_x[i,j,:] = e0_1d
eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)


    

    




  File "<ipython-input-34-f14ef29654fd>", line 12
    eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)
                                                       ^
SyntaxError: invalid syntax

In [ ]: