In [30]:

    

from SimPEG import *
import simpegEM as EM
from simpegem1d import *
from scipy.constants import mu_0
%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [239]:

    

meshType = 'CYL'
cs, ncx, ncz, npad = 10., 10, 5, 15
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')


    

In [240]:

    

active = mesh.vectorCCz<0.
layer1 = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-60.)
layer2 = (mesh.vectorCCz<-60) & (mesh.vectorCCz>=-100.)
layer3 = (mesh.vectorCCz<-100) & (mesh.vectorCCz>=-200.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap


    

In [241]:

    

sig_half = 1e-3
sigma = np.ones(mesh.nCz)*sig_air
sigma[active] = sig_half
mtrue = np.log(sigma[active])


    

In [242]:

    

frequency = np.logspace(2, 5, 31)
xyz_rx = np.r_[8., 0., 30.]
rx1 = EM.FDEM.RxFDEM(xyz_rx.reshape([1,-1]), 'bzr') 
rx2 = EM.FDEM.RxFDEM(xyz_rx.reshape([1,-1]), 'bzi')
txList = []
for freq in frequency:    
    tx = EM.FDEM.TxFDEM(np.array([0., 0., 30.]), 'VMD', freq, [rx1, rx2])
    tx.radius = 100.
    txList.append(tx)
    
survey = EM.FDEM.SurveyFDEM(txList)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=False)
survey.pair(prb)


    

In [243]:

    

Bz_simpeg = survey.dpred(mtrue)


    

In [244]:

    

FDsurvey = BaseEM1D.EM1DSurveyFD()
FDsurvey.rxLoc = np.array([0., 0., 30.])
FDsurvey.txLoc = np.array([0., 0., 30.])
FDsurvey.fieldtype = 'secondary'
nearthick = np.logspace(-1, 1, 5)
deepthick = np.logspace(1, 2, 10)
hx = np.r_[nearthick, deepthick]
mesh1D = Mesh.TensorMesh([hx], [0.])
depth = -mesh1D.gridN[:-1]
LocSigZ = -mesh1D.gridCC
nlay = depth.size
topo = np.r_[0., 0., 0.]
FDsurvey.depth = depth
FDsurvey.topo = topo
FDsurvey.LocSigZ = LocSigZ
FDsurvey.frequency = frequency
FDsurvey.Nfreq = FDsurvey.frequency.size
FDsurvey.Setup1Dsystem()        
expmap = BaseEM1D.BaseEM1DMap(mesh1D)
modelReal = expmap
m_1D = np.log(np.ones(nlay)*sig_half)
FDsurvey.rxType = 'Hz'
FDsurvey.switchRI = 'all'
WT0, WT1, YBASE = DigFilter.LoadWeights()
options = {'WT0': WT0, 'WT1': WT1, 'YBASE': YBASE}
prob = EM1D.EM1D(mesh1D, modelReal, **options)
prob.pair(FDsurvey)
prob.chi = np.zeros(FDsurvey.nlay)


    

In [245]:

    

prob.CondType = 'Real'
# prob.survey.txType = 'CircularLoop'
# I = 1e0
# a = 1e2
# prob.survey.I = I
# prob.survey.a = a
prob.survey.txType = 'VMD'
prob.survey.offset = 8
Bz_anal = FDsurvey.dpred(m_1D)


    

In [246]:

    

print Hz_anal


    

    




[ -1.24388034e-16  -1.56040332e-16  -1.95703951e-16  -2.45392432e-16
  -3.07620759e-16  -3.85528733e-16  -4.83034211e-16  -6.05023282e-16
  -7.57586085e-16  -9.48308943e-16  -1.18663592e-15  -1.48431579e-15
  -1.85595410e-15  -2.31969401e-15  -2.89805518e-15  -3.61896579e-15
  -4.51703047e-15  -5.63508571e-15  -7.02610490e-15  -8.75552746e-15
  -1.09041016e-14  -1.35713468e-14  -1.68797632e-14  -2.09799362e-14
  -2.60567142e-14  -3.23366642e-14  -4.00970462e-14  -4.96765844e-14
  -6.14883542e-14  -7.60351496e-14  -9.39277414e-14  -1.01343469e-14
  -1.18037768e-14  -1.37471027e-14  -1.60089766e-14  -1.86412566e-14
  -2.17041513e-14  -2.52675410e-14  -2.94125017e-14  -3.42330621e-14
  -3.98382266e-14  -4.63543023e-14  -5.39275735e-14  -6.27273706e-14
  -7.29495873e-14  -8.48207063e-14  -9.86023986e-14  -1.14596770e-13
  -1.33152334e-13  -1.54670796e-13  -1.79614744e-13  -2.08516346e-13
  -2.41987146e-13  -2.80729087e-13  -3.25546844e-13  -3.77361596e-13
  -4.37226313e-13  -5.06342672e-13  -5.86079643e-13  -6.77993806e-13
  -7.83851396e-13  -9.05652020e-13]

In [247]:

    

Bz_Simpeg = Bz_simpeg.reshape((2,FDsurvey.Nfreq), order='F')
plt.loglog(frequency, -Bz_anal[:FDsurvey.Nfreq], 'k-')
plt.loglog(frequency, -Bz_anal[FDsurvey.Nfreq:], 'b-')
plt.loglog(FDsurvey.frequency, -Bz_Simpeg[0,:], 'ko')
plt.loglog(FDsurvey.frequency, -Bz_Simpeg[1,:], 'bo')


    

    
Out[247]:





[<matplotlib.lines.Line2D at 0xff19a50>]