Numerical examples of the book
Anton Ziolkowski and Evert Slob, 2019, Cambridge University Press; ISBN: 9781107058620.
Copyright 2018 Dieter Werthmüller
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
import numpy as np
import empymod as epm
import matplotlib.pyplot as plt
%matplotlib inline
The book shows in the chapter black-and-white figures, and in the plate section coloured versions.
plt.rcParams['image.cmap'] = 'viridis' # Colour
# plt.rcParams['image.cmap'] = 'gray' # Black-and-white
# Offsets for period-amplitude plots
x = 5e2+np.linspace(0, 1.95e4, 196)
# Periods for offset-amplitude plots
p = np.logspace(-1.5, 1, 401)
# Common model parameters
model = {'src': [0, 0, 1e-3], # Source location
'rec': [x, x*0, 1e-3], # Receiver location
'freqtime': 1/p, # Frequencies
'res': [2e14, 20, 400, 20], # Resistivity
'aniso': [1, np.sqrt(2), np.sqrt(2), np.sqrt(2)], # Anisotropy
'opt': 'parallel', # Optimisation
'loop': 'off', # Loop over frequencies
'htarg': 'key_201_2012', # Use same filter as in book
'verb': 1} # Verbosity
# Target at 1 km depth
ex1km = epm.dipole(depth=[0, 1000, 1040], **model)
# Target at 3 km depth
ex3km = epm.dipole(depth=[0, 3000, 3040], **model)
# Analytical, diffusive half-space solution
hs = epm.analytical(src=[0, 0, 0], rec=[x, x*0, 0], res=20, aniso=np.sqrt(2),
freqtime=1/p, solution='dhs', verb=1)
def plot_result(data, error=False):
"""Create figure.
Most figures in this section are very similar, and just differ in
the shown data. We can therefore define a plot-function which
we subsequently call with the different results.
"""
n = len(data)
# Define figure size depending if 2 or 4 data-sets.
if n == 2:
tit = 'Depth to resistor: '
titadd = ['1 km', '3 km']
size = (12, 6)
clbargs = {'location': 'bottom', 'fraction': .05,
'pad': 0.2, 'aspect': 30}
else:
tit = 'Half-space response'
titadd = ['', '']
size = (6, 6)
clbargs = {'location': 'bottom'}
if error:
clbticks = [-2, -1.5, -1, -0.5, 0, 0.5, 1]
vmin = -2
vmax = 1
else:
clbticks = [-13, -13, -12, -11, -10, -9]
vmin = -13
vmax = -9
# Start figure
fig, axs = plt.subplots(figsize=size, nrows=1, ncols=n)
if n == 2:
axs = axs.ravel()
plt.subplots_adjust(hspace=0.3, wspace=0.3)
else:
axs = [axs,]
# Loop over data
for i, val in enumerate(data):
plt.sca(axs[i])
plt.title(tit + titadd[i % 2], fontweight='bold')
# pcolormesh takes x/y as start and end point of pixels,
# not midpoints. So we have to create these.
dx = (x[1]-x[0])/2
px = np.r_[x-dx, x[-1]+dx]
lp = np.log10(p)
dp = (lp[1]-lp[0])/2
pp = 10**(np.r_[lp-dp, lp[-1]+dp])
# Plot result
cs = plt.pcolormesh(px/1000, pp, np.log10(np.abs(val)),
vmin=vmin, vmax=vmax)
# Plot contours
level_factor = 2
if (vmax-vmin) < 4:
level_factor *= 2
levels = np.linspace(vmin, vmax, level_factor*(vmax-vmin)+1)
cs2 = plt.contour(x/1000, p, np.log10(np.abs(val)),
levels=levels,
linewidths=0.5, colors=[(0, 0, 0, 0.5)])
plt.xticks([5, 10, 15, 20])
plt.xlabel('Offset (km)')
plt.ylabel(r'$\log_{10}$(period) (s)')
axs[i].set_yscale('log')
axs[i].invert_yaxis()
plt.ylim([np.max(p), np.min(p)])
plt.yticks([1e-1, 1e0, 1e1], ('-1', '0', '1'))
# Plot colorbar
cax, kw = plt.matplotlib.colorbar.make_axes(axs, **clbargs)
cb = plt.colorbar(cs, cax=cax, ticks=clbticks, **kw)
plt.show()
def set_axis(ytype, p=False):
"""General settings for line plots."""
if np.any(p):
plt.xlabel('Period (s)')
if ytype == 1:
plt.ylabel('Electric field (pV/m)')
else:
plt.ylabel('Phase (rad)')
plt.xlim([np.min(p), np.max(p)])
else:
plt.xlabel('Offset (km)')
if ytype == 1:
plt.ylabel('Electric field (V/m)')
else:
plt.ylabel('Phase (rad)')
plt.xlim([0, 20])
plt.xticks([0, 5, 10, 15, 20])
plt.legend()
def thin_y_ticks():
"""Hide every second y-tick label."""
for label in plt.gca().get_yticklabels()[::2]:
label.set_visible(False)
plot_result([hs])
Figure 5.46 Electric field amplitude for the VTI half-space land model as a function of offset and source oscillation period. The colour bar indicates the values of the logarithm of the amplitude.
plot_result([ex1km, ex3km])
Figure 5.47 Electric field amplitudes for the land VTI model with the buried resistive layer at 1 km depth (left) and at 3 km depth (right) as a function of offset and source oscillation period. The colour bar indicates the amplitudes on a logarithmic scale.
plot_result([(ex1km-hs)/hs, (ex3km-hs)/hs], True)
Figure 5.48 Normalised amplitude of the difference in the electric fields for the land VTI model with a buried resistive layer at 1 km depth (left) and at 3 km depth (right) as a function of offset and source oscillation period. The colour bar indicates the amplitude of the normalised difference in the electric field on a logarithmic scale.
In red is the analytical half-space solution, in black the solution from the modeller.
# Define indices
nx1, nx3 = 55, 158
np1, np3 = 177, 317
plt.figure(figsize=(10, 8))
plt.subplots_adjust(wspace=.3, hspace=.4)
# 1st subplot
plt.subplot(221)
plt.title('Period = '+str(np.round(p[np1], 1))+' s; depth to resistor: 1 km',
fontweight='bold')
plt.semilogy(x/1000, np.abs(ex1km[np1, :]), 'k-', label='With resistor')
plt.semilogy(x/1000, np.abs(hs[np1, :]), 'k--', label='Without resistor')
set_axis(1)
thin_y_ticks()
# 2nd subplot
plt.subplot(223)
plt.title('Period = '+str(np.round(p[np3], 1))+' s; depth to resistor: 3 km',
fontweight='bold')
plt.semilogy(x/1000, np.abs(ex3km[np3, :]), 'k-', label='With resistor')
plt.semilogy(x/1000, np.abs(hs[np3, :]), 'k--', label='Without resistor')
set_axis(1)
thin_y_ticks()
# 3rd subplot
plt.subplot(222)
plt.title('Offset = '+str(np.round(x[nx1])/1000) +
' km; depth to resistor: 1 km', fontweight='bold')
plt.semilogx(p, np.abs(ex1km[:, nx1])/1e-12, 'k-', label='With resistor')
plt.semilogx(p, np.abs(hs[:, nx1])/1e-12, 'k--', label='Without resistor')
set_axis(1, p)
# 4th subplot
plt.subplot(224)
plt.title('Offset = '+str(np.round(x[nx3])/1000) +
' km; depth to resistor: 3 km', fontweight='bold')
plt.semilogx(p, np.abs(ex3km[:, nx3])/1e-12, 'k-', label='With resistor')
plt.semilogx(p, np.abs(hs[:, nx3])/1e-12, 'k--', label='Without resistor')
set_axis(1, p)
plt.show()
Figure 5.49 Electric field amplitudes for the land VTI model with (solid lines) and without (dashed lines) a buried resistive layer at 1 km depth (top) and at 3 km depth (bottom) as a function of offset (left) and source oscillation period (right).
plt.figure(figsize=(10, 8))
plt.subplots_adjust(wspace=.3, hspace=.4)
# 1st subplot
plt.subplot(221)
plt.title('Period = '+str(np.round(p[np1], 1))+' s; depth to resistor: 1 km',
fontweight='bold')
plt.plot(x/1000, -np.angle(ex1km[np1, :]), 'k-', label='With resistor')
plt.plot(x/1000, -np.angle(hs[np1, :]), 'k--', label='Without resistor')
set_axis(2)
# 2nd subplot
plt.subplot(223)
plt.title('Period = '+str(np.round(p[np3], 1))+' s; depth to resistor: 3 km',
fontweight='bold')
plt.plot(x/1000, -np.angle(ex3km[np3, :]), 'k-', label='With resistor')
plt.plot(x/1000, -np.angle(hs[np3, :]), 'k--', label='Without resistor')
set_axis(2)
# 3rd subplot
plt.subplot(222)
plt.title('Offset = '+str(np.round(x[nx1])/1000) +
' km; depth to resistor: 1 km', fontweight='bold')
plt.semilogx(p, -np.angle(ex1km[:, nx1]), 'k-', label='With resistor')
plt.semilogx(p, -np.angle(hs[:, nx1]), 'k--', label='Without resistor')
set_axis(2, p)
# 4th subplot
plt.subplot(224)
plt.title('Offset = '+str(np.round(x[nx3])/1000) +
' km; depth to resistor: 3 km', fontweight='bold')
plt.semilogx(p, -np.angle(ex3km[:, nx3]), 'k-', label='With resistor')
plt.semilogx(p, -np.angle(hs[:, nx3]), 'k--', label='Without resistor')
set_axis(2, p)
plt.show()
Figure 5.50 Electric field phases for the land VTI model with (solid lines) and without (dashed lines) a buried resistive layer at 1 km depth (top) and at 3 km depth (bottom) as a function of offset (left) and source oscillation period (right).
epm.versions('HTML')