In [1]:

    

from SimPEG import *
import simpegEM as EM
from simpegem1d import Utils1D
%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 100


    

In [3]:

    

matplotlib.rcParams.update({'font.size': 16})


    

In [4]:

    

meshType = 'CYL'
cs, ncx, ncz, npad = 20., 25, 30, 12
hx = [(cs,ncx), (cs,npad,1.3)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.CylMesh([hx,1,hz], '00C')


    

In [5]:

    

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


    

In [6]:

    

sig_half = 1e-3
sig_air = 1e-8
sig_layer1 = 1./300
sig_layer2 = 1./5.
sigma = np.ones(mesh.nCz)*sig_air
sigma_half = sigma.copy()
sigma[active] = sig_half
sigma[layer1] = sig_layer1
sigma[layer2] = sig_layer2
sigma_half[active] = sig_half
sigma_half[layer1] = sig_layer1
sigma_half[layer2] = 1./500.
mtrue = np.log(sigma[active])


    

In [7]:

    

sig_sea = (sigma[active])
sig_half = sigma_half[active]


    

In [8]:

    

fig, ax = plt.subplots(1,1, figsize = (6, 6))
ax.plot(1000, 1000, 'k')
ax.plot(1000, 1000, 'b')
ax.legend(('With intrusion', 'Without intrusion'), loc=4)
Utils1D.plotLayer(sig_half, mesh.vectorCCz[active], ax = ax, **{'color':'b'})
Utils1D.plotLayer(sig_sea, mesh.vectorCCz[active], showlayers=True, ax = ax)
ax.set_ylim(-500., 0.)
fig.savefig('./figures/1Dmodel.png', dpi=200)


    

In [9]:

    

prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=False)


    

In [10]:

    

prb.Solver = SolverLU


    

In [11]:

    

frequency = np.logspace(-1, 3, 41)


    

In [12]:

    

x = np.linspace(0, 250, 26)
xyz_rx = Utils.ndgrid(x, np.r_[0.], np.r_[0.])


    

In [13]:

    

rxr = EM.FDEM.RxFDEM(xyz_rx,'bzr')
rxi = EM.FDEM.RxFDEM(xyz_rx,'bzi')
txList = []
for freq in frequency:
    tx = EM.FDEM.TxFDEM(np.array([0., 0., 0.]), 'CircularLoop', freq, [rxr, rxi])
    tx.radius = 250.
    txList.append(tx)

survey = EM.FDEM.SurveyFDEM(txList)
if prb.ispaired:
    prb.unpair()
if survey.ispaired:
    survey.unpair()
prb.pair(survey)


    

In [14]:

    

bz_sea = survey.dpred(np.log(sig_sea))


    

In [15]:

    

bz_half = survey.dpred(np.log(sig_half))


    

In [16]:

    

Bz_sea = bz_sea.reshape((x.size, 2, frequency.size), order='F')
Bz_half = bz_half.reshape((x.size, 2, frequency.size), order='F')


    

In [17]:

    

Bz_sea_complex = (Bz_sea[:,0,:]+1j*Bz_sea[:,1,:]).reshape((x.size, frequency.size), order='F')
Bz_half_complex = (Bz_half[:,0,:]+1j*Bz_half[:,1,:]).reshape((x.size, frequency.size), order='F')


    

In [18]:

    

dum = np.log(Bz_sea_complex/Bz_half_complex)


    

In [19]:

    

Bz_ratio = np.exp(dum.real)
Bz_phasedif = dum.imag


    

In [20]:

    

FreqX = Utils.ndgrid(x, frequency)
X = FreqX[:,0].reshape((x.size, frequency.size), order='F')
Freq = FreqX[:,1].reshape((x.size, frequency.size), order='F')


    

In [21]:

    

print ('$10^{%3.2f}$') % (3.0)


    

    




$10^{3.00}$

In [22]:

    

test = (np.arange(9)*0.5-2).tolist()
dummy = []
for itest in test:
    dummy.append('10$^{'+str(itest)+'}$')


    

In [24]:

    

fig, ax = plt.subplots(1,1, figsize = (8,6))
vmin = Bz_ratio.min()
vmax = Bz_ratio.max()
dat = ax.contourf(X, np.log10(Freq), (Bz_ratio), 100, vmin=vmin, vmax=vmax, clim=(vmin, vmax), cmap = cm.RdPu)
ax.set_xlabel('Distance (m)', fontsize = 16)
ax.set_ylabel('Frequency (Hz)', fontsize = 16)
ax.set_yticklabels(dummy)
cb = colorbar(dat, ax = ax, format='%5.2f', ticks=np.linspace(vmin, vmax, 5))
levels = np.r_[5., 10., 20.]
CS = ax.contour(X, np.log10(Freq), (Bz_ratio), levels, colors='k')
ax.clabel(CS, inline=1, fmt = '%3.1f')
cb.set_label('Amplitude ratio', fontsize = 16)
fig.savefig('figures/ampratio_freq.png', dpi=200)


    

In [25]:

    

fig, ax = plt.subplots(1,1, figsize = (8,6))
vmin = Bz_phasedif.min()*1e3
vmax = Bz_phasedif.max()*1e3
dat = ax.contourf(X, np.log10(Freq), 1e3*(Bz_phasedif), 100, vmin=vmin, vmax=vmax, clim=(vmin, vmax), colormaps=cm.jet_r)
ax.set_xlabel('Distance (m)', fontsize = 16)
ax.set_ylabel('Frequency (Hz)', fontsize = 16)
ax.set_yticklabels(dummy)
cb = colorbar(dat, ax = ax, format='%5.1f', ticks=np.linspace(vmin, vmax, 5))
cb.set_label('Phase difference (mrad)', fontsize = 14)
fig.savefig('figures/phasedif_freq.png', dpi=200)