"""
Part of comet/pyhed/hed
"""
# COMET - COupled Magnetic resonance Electrical resistivity Tomography
# Copyright (C) 2019 Nico Skibbe
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from __future__ import print_function
import numpy as np
from comet.pyhed.misc import vec
from comet.pyhed import log
# import time
# import sys
from . libHED import World1D, HED
# from comet import pyhed as ph
[docs]def hankelfc(order):
""" Filter coefficients for hankel transformation by Anderson (1980) """
if order == 1: # sin
fc = np.array([
2.59526236e-7, 3.66544843e-7, 5.17830795e-7, 7.31340622e-7,
1.03322805e-6, 1.45918500e-6, 2.06161065e-6, 2.91137793e-6,
4.11357863e-6, 5.80876420e-6, 8.20798075e-6, 1.15895083e-5,
1.63778560e-5, 2.31228459e-5, 3.26800649e-5, 4.61329334e-5,
6.52101085e-5, 9.20390575e-5, 1.30122935e-4, 1.83620431e-4,
2.59656626e-4, 3.66311982e-4, 5.18141184e-4, 7.30717340e-4,
1.03392184e-3, 1.45742714e-3, 2.06292302e-3, 2.90599911e-3,
4.11471902e-3, 5.79042763e-3, 8.20004722e-3, 1.15192930e-2,
1.63039133e-2, 2.28257757e-2, 3.22249222e-2, 4.47864328e-2,
6.27329625e-2, 8.57059100e-2, 1.17418314e-1, 1.53632655e-1,
1.97717964e-1, 2.28849849e-1, 2.40311038e-1, 1.65409220e-1,
2.84701476e-3, -2.88016057e-1, -3.69097406e-1, -2.50107514e-2,
5.71811256e-1, -3.92261572e-1, 7.63280044e-2, 5.16233994e-2,
-6.48012082e-2, 4.89047141e-2, -3.26936331e-2, 2.10539842e-2,
-1.33862549e-2, 8.47124695e-3, -5.35123972e-3, 3.37796651e-3,
-2.13174466e-3, 1.34513833e-3, -8.48749612e-4, 5.35531006e-4,
-3.37898780e-4, 2.13200109e-4, -1.34520273e-4, 8.48765787e-5,
-5.35535069e-5, 3.37899801e-5, -2.13200365e-5, 1.34520337e-5,
-8.48765949e-6, 5.35535110e-6, -3.37899811e-6, 2.13200368e-6,
-1.34520338e-6, 8.48765951e-7, -5.35535110e-7, 3.37899811e-7],
np.float)
nc = np.int(80)
nc0 = np.int(40)
elif order == 2: # cos
fc = np.array([
1.63740363e-7, 1.83719709e-7, 2.06136904e-7, 2.31289411e-7,
2.59510987e-7, 2.91176117e-7, 3.26704977e-7, 3.66569013e-7,
4.11297197e-7, 4.61483045e-7, 5.17792493e-7, 5.80972733e-7,
6.51862128e-7, 7.31401337e-7, 8.20645798e-7, 9.20779729e-7,
1.03313185e-6, 1.15919300e-6, 1.30063594e-6, 1.45933752e-6,
1.63740363e-6, 1.83719709e-6, 2.06136904e-6, 2.31289411e-6,
2.59510987e-6, 2.91176117e-6, 3.26704977e-6, 3.66569013e-6,
4.11297197e-6, 4.61483045e-6, 5.17792493e-6, 5.80972733e-6,
6.51862128e-6, 7.31401337e-6, 8.20645798e-6, 9.20779729e-6,
1.03313185e-5, 1.15919300e-5, 1.30063594e-5, 1.45933752e-5,
1.63740363e-5, 1.83719709e-5, 2.06136904e-5, 2.31289411e-5,
2.59510987e-5, 2.91176117e-5, 3.26704977e-5, 3.66569013e-5,
4.11297197e-5, 4.61483045e-5, 5.17792493e-5, 5.80972733e-5,
6.51862128e-5, 7.31401337e-5, 8.20645798e-5, 9.20779729e-5,
1.03313185e-4, 1.15919300e-4, 1.30063594e-4, 1.45933752e-4,
1.63740363e-4, 1.83719709e-4, 2.06136904e-4, 2.31289411e-4,
2.59510987e-4, 2.91176117e-4, 3.26704976e-4, 3.66569013e-4,
4.11297197e-4, 4.61483045e-4, 5.17792493e-4, 5.80972733e-4,
6.51862127e-4, 7.31401337e-4, 8.20645797e-4, 9.20779730e-4,
1.03313185e-3, 1.15919300e-3, 1.30063593e-3, 1.45933753e-3,
1.63740362e-3, 1.83719710e-3, 2.06136901e-3, 2.31289411e-3,
2.59510977e-3, 2.91176115e-3, 3.26704948e-3, 3.66569003e-3,
4.11297114e-3, 4.61483003e-3, 5.17792252e-3, 5.80972566e-3,
6.51861416e-3, 7.31400728e-3, 8.20643673e-3, 9.20777603e-3,
1.03312545e-2, 1.15918577e-2, 1.30061650e-2, 1.45931339e-2,
1.63734419e-2, 1.83711757e-2, 2.06118614e-2, 2.31263461e-2,
2.59454421e-2, 2.91092045e-2, 3.26529302e-2, 3.66298115e-2,
4.10749753e-2, 4.60613861e-2, 5.16081994e-2, 5.78193646e-2,
6.46507780e-2, 7.22544422e-2, 8.03873578e-2, 8.92661837e-2,
9.80670729e-2, 1.07049506e-1, 1.13757572e-1, 1.18327217e-1,
1.13965041e-1, 1.00497783e-1, 6.12958082e-2, -1.61234222e-4,
-1.11788551e-1, -2.27536948e-1, -3.39004453e-1, -2.25128800e-1,
8.98279919e-2, 5.12510388e-1, -1.31991937e-1, -3.35136479e-1,
3.64868100e-1, -2.34039961e-1, 1.32085237e-1, -7.56739672e-2,
4.52296662e-2, -2.78297002e-2, 1.73727753e-2, -1.09136894e-2,
6.87397283e-3, -4.33413470e-3, 2.73388730e-3, -1.72477355e-3,
1.08821012e-3, -6.86602007e-4, 4.33213523e-4, -2.73338487e-4,
1.72464733e-4, -1.08817842e-4, 6.86594042e-5, -4.33211523e-5,
2.73337984e-5, -1.72464607e-5, 1.08817810e-5, -6.86593962e-6,
4.33211503e-6, -2.73337979e-6, 1.72464606e-6, -1.08817810e-6,
6.86593961e-7, -4.33211503e-7, 2.73337979e-7, -1.72464606e-7],
np.float)
nc = np.int(164)
nc0 = np.int(122)
elif order == 3: # J0
fc = np.array([
2.89878288e-7, 3.64935144e-7, 4.59426126e-7, 5.78383226e-7,
7.28141338e-7, 9.16675639e-7, 1.15402625e-6, 1.45283298e-6,
1.82900834e-6, 2.30258511e-6, 2.89878286e-6, 3.64935148e-6,
4.59426119e-6, 5.78383236e-6, 7.28141322e-6, 9.16675664e-6,
1.15402621e-5, 1.45283305e-5, 1.82900824e-5, 2.30258527e-5,
2.89878259e-5, 3.64935186e-5, 4.59426051e-5, 5.78383329e-5,
7.28141144e-5, 9.16675882e-5, 1.15402573e-4, 1.45283354e-4,
1.82900694e-4, 2.30258630e-4, 2.89877891e-4, 3.64935362e-4,
4.59424960e-4, 5.78383437e-4, 7.28137738e-4, 9.16674828e-4,
1.15401453e-3, 1.45282561e-3, 1.82896826e-3, 2.30254535e-3,
2.89863979e-3, 3.64916703e-3, 4.59373308e-3, 5.78303238e-3,
7.27941497e-3, 9.16340705e-3, 1.15325691e-2, 1.45145832e-2,
1.82601199e-2, 2.29701042e-2, 2.88702619e-2, 3.62691810e-2,
4.54794031e-2, 5.69408192e-2, 7.09873072e-2, 8.80995426e-2,
1.08223889e-1, 1.31250483e-1, 1.55055715e-1, 1.76371506e-1,
1.85627738e-1, 1.69778044e-1, 1.03405245e-1, -3.02583233e-2,
-2.27574393e-1, -3.62173217e-1, -2.05500446e-1, 3.37394873e-1,
3.17689897e-1, -5.13762160e-1, 3.09130264e-1, -1.26757592e-1,
4.61967890e-2, -1.80968674e-2, 8.35426050e-3, -4.47368304e-3,
2.61974783e-3, -1.60171357e-3, 9.97717882e-4, -6.26275815e-4,
3.94338818e-4, -2.48606354e-4, 1.56808604e-4, -9.89266288e-5,
6.24152398e-5, -3.93805393e-5, 2.48472358e-5, -1.56774945e-5,
9.89181741e-6, -6.24131160e-6, 3.93800058e-6, -2.48471018e-6,
1.56774609e-6, -9.89180896e-7, 6.24130948e-7, -3.93800005e-7,
2.48471005e-7, -1.56774605e-7, 9.89180888e-8, -6.24130946e-8],
np.float)
nc = np.int(100)
nc0 = np.int(60)
elif order == 4: # J1
fc = np.array([
1.84909557e-13, 2.85321327e-13, 4.64471808e-13, 7.16694771e-13,
1.16670043e-12, 1.80025587e-12, 2.93061898e-12, 4.52203829e-12,
7.36138206e-12, 1.13588466e-11, 1.84909557e-11, 2.85321327e-11,
4.64471808e-11, 7.16694771e-11, 1.16670043e-10, 1.80025587e-10,
2.93061898e-10, 4.52203829e-10, 7.36138206e-10, 1.13588466e-9,
1.84909557e-9, 2.85321326e-9, 4.64471806e-9, 7.16694765e-9,
1.16670042e-8, 1.80025583e-8, 2.93061889e-8, 4.52203807e-8,
7.36138149e-8, 1.13588452e-7, 1.84909521e-7, 2.85321237e-7,
4.64471580e-7, 7.16694198e-7, 1.16669899e-6, 1.80025226e-6,
2.93060990e-6, 4.52201549e-6, 7.36132477e-6, 1.13587027e-5,
1.84905942e-5, 2.85312247e-5, 4.64449000e-5, 7.16637480e-5,
1.16655653e-4, 1.79989440e-4, 2.92971106e-4, 4.51975783e-4,
7.35565435e-4, 1.13444615e-3, 1.84548306e-3, 2.84414257e-3,
4.62194743e-3, 7.10980590e-3, 1.15236911e-2, 1.76434485e-2,
2.84076233e-2, 4.29770596e-2, 6.80332569e-2, 9.97845929e-2,
1.51070544e-1, 2.03540581e-1, 2.71235377e-1, 2.76073871e-1,
2.16691977e-1, -7.83723737e-2, -3.40675627e-1, -3.60693673e-1,
5.13024526e-1, -5.94724729e-2, -1.95117123e-1, 1.99235600e-1,
-1.38521553e-1, 8.79320859e-2, -5.50697146e-2, 3.45637848e-2,
-2.17527180e-2, 1.37100291e-2, -8.64656417e-3, 5.45462758e-3,
-3.44138864e-3, 2.17130686e-3, -1.36998628e-3, 8.64398952e-4,
-5.45397874e-4, 3.44122545e-4, -2.17126585e-4, 1.36997597e-4,
-8.64396364e-5, 5.45397224e-5, -3.44122382e-5, 2.17126544e-5,
-1.36997587e-5, 8.64396338e-6, -5.45397218e-6, 3.44122380e-6,
-2.17126543e-6, 1.36997587e-6, -8.64396337e-7, 5.45397218e-7],
np.float)
nc = np.int(100)
nc0 = np.int(60)
return (np.reshape(fc, (-1, 1)), nc, nc0) # (100,) -> (100, 1)
[docs]def downward(u, model, rho, d, f, mode):
""" Downward attenuation.
Parameters
----------
u: np.ndarray
Horizontal wavenumbers based on Hankel factors and horizontal
tx-rx distance. Shape: (Hankel, n_rx)
model: np.ndarray
Polar coords of the receiver pos (3, n).
rho: np.ndarray
Resistivities for each layer (Ohm*m).
d: np.ndarray
Thicknesses of each layer (m).
f: float
Frequency (Hz).
Returns
-------
aa : np.ndarray
Ratio of the partial wave amplitude A(z,u)/A(0,u)
aap : np.andarray
Ratio of the partial wave amplitude A'(z,u)/A'(0,u)
bt : np.ndarray
Admittance at the surface of the layerd halfspace
"""
mu0 = 4e-7 * np.pi
iwm = 1j * 2. * np.pi * f * mu0
nl = rho.size
h = np.zeros(nl)
h[1:] += np.cumsum(d)
rho = rho[np.newaxis, :]
rho = np.tile(rho.T, (1, u.shape[0]))
rho = rho[:, :, np.newaxis]
d = d[np.newaxis, :]
d = np.tile(d.T, (1, u.shape[0]))
d = d[:, :, np.newaxis]
alpha = np.sqrt(u ** 2 + iwm/rho) # dim: Rho x Hankel x model
if nl == 1: # homogeneous halfspace
b1 = alpha[0]
a = np.exp((-alpha[0]) * (-model[2]))
ap = a.copy()
return b1, a, ap
elif nl > 1: # multi-layer case
ealphad = np.exp(-2.0 * alpha[:-1] * d) # needs time
talphad = (1.0 - ealphad) / (ealphad + 1.0)
# print(min(-alpha[:-1] * d), max(-alpha[:-1] * d))
# save for later
b = np.copy(alpha)
# recursive calculation from bottom to top for nl-1 layer
beta = np.ones_like(alpha)
if mode.lower() == 'tm':
beta[:-1] *= rho[:-1]/rho[1:]
albe = alpha[:-1] * beta[:-1]
elif mode.lower() == 'te':
albe = alpha[:-1]
else:
raise Exception('{} not a valid mode for calculation of reflection\
coefficients. Choose either "te" or "tm".'.format(mode))
albetal = albe * talphad
for n in range(nl-2, -1, -1):
b[n] *= (b[n + 1] + albetal[n]) / (
albe[n] + b[n + 1] * talphad[n])
# Impedance
c = 1.0 / b
b1 = np.copy(b[0]) # only for surface
# Variation from one layer boundary to the other
dim2 = u.shape[0]
dim3 = u.shape[1]
aa = np.zeros((nl-1) * dim2 * dim3, dtype=np.complex).reshape(
nl-1, dim2, dim3)
aap = np.zeros_like(aa)
for n in range(0, nl-1):
# note: here comes the memory
term_aa = beta[n] * (b[n] + alpha[n]) / (b[n+1] + albe[n])
expterm = np.exp(-alpha[n] * d[n])
# aa[n] = beta[n] * (b[n] + alpha[n]) / \
# (b[n+1] + albe[n]) * \
# np.exp(-alpha[n] * d[n])
term_aap = (alpha[n] * c[n] + 1.0) / (albe[n] * c[n+1] + 1.0)
# aap[n] = (alpha[n] * c[n] + 1.0) / (1.0 + albe[n] *
# c[n+1]) * \
# np.exp(-alpha[n] * d[n])
aa[n] = term_aa * expterm
aap[n] = term_aap * expterm
# del bla
# del bla2
# del expterm
# find corresponding layer for each node/receiver
ind = np.zeros(np.shape(model)[1], dtype=np.int)
for i in range(0, len(h) - 1):
ind[np.nonzero(np.logical_and(model[2] <= -h[i],
model[2] > -h[i + 1]))] += i
ind[np.nonzero(model[2] <= -h[-1])] = len(h) - 1
a = np.zeros(np.shape(u), dtype=np.complex128)
ap = np.zeros(np.shape(u), dtype=np.complex128)
for i in range(0, nl):
# all layer at once (except the last)
idl = np.where(ind == i)[0]
if i < (nl - 1):
prod_aa = np.prod(aa[0:i, :, idl], 0)
prod_aap = np.prod(aap[0:i, :, idl], 0)
# aap[1:n-1] = sum[1:n-1](F'(n + 1) / F'(n))
# aa[1:n-1] = sum[1:n-1](F(n + 1) / F(n))
exp_1 = np.exp(-alpha[i, :, idl].T * (-model[2][idl] - h[i]))
exp_2 = np.exp(-alpha[i, :, idl].T *
(d[i, :, :] + h[i+1] + model[2][idl]))
factor_a = 0.5 * (
1.0 + b[i, :, idl].T / alpha[i, :, idl].T) * (
exp_1 -
(b[i+1, :, idl].T -
albe[i, :, idl].T) /
(b[i+1, :, idl].T +
albe[i, :, idl].T) *
exp_2)
factor_ap = 0.5 * \
(1.0 + alpha[i, :, idl].T * c[i, :, idl].T) * \
(exp_1 +
(1.0 - albe[i, :, idl].T * c[i+1, :, idl].T) /
(1.0 + albe[i, :, idl].T * c[i+1, :, idl].T) *
exp_2)
a[:, idl] = prod_aa * factor_a
ap[:, idl] = prod_aap * factor_ap
else:
# last layer (lower halfspace)
exp_alphai = np.exp(-alpha[i, :, idl].T *
(-model[2][idl] - h[i]))
a[:, idl] = np.prod(aa[:, :, idl], 0) * exp_alphai
ap[:, idl] = np.prod(aap[:, :, idl], 0) * exp_alphai
return b1, a, ap
[docs]def btp(u, model, rho, d, f, mode):
""" Airspace only, internal function, for imput see **downout**.
Do not call directly.
"""
nl = len(rho) # number of layer
mu0 = 4e-7 * np.pi # mu null
iwm = 1j * 2. * np.pi * f * mu0
uiwrho = np.sqrt((u**2)[:, :, np.newaxis] + iwm/rho[:nl])
b = uiwrho[:, :, -1]
alpha = uiwrho[:, :, :-1]
if mode == 'te':
beta = np.ones(nl - 1)
elif mode == 'tm':
beta = rho[:nl-1]/rho[1:nl]
else:
raise Exception('"{}" not a valid mode for calculation of reflection\
coefficients. Choose either "te" or "tm".'.format(mode))
cth_exp = np.exp(-2. * alpha * d[:nl-1])
cth = (1.0 - cth_exp) / (cth_exp + 1.0)
# cth= np.tanh(-alpha * d[:nl-1]) # note: slower, less accurate
if nl > 1:
b_last = b
for i in range(nl - 2, -1, -1): # second-to-last layer to top layer
b_new = (b_last + alpha[:, :, i] * beta[i] * cth[:, :, i]) / \
(beta[i] + cth[:, :, i] * b_last / alpha[:, :, i])
b_last = b_new
return b_new
else:
return b
[docs]def downout(u, model, rho, d, f, mode):
""" Overall call function for recursive calculation.
Parameters
----------
u: np.ndarray
Horizontal wavenumbers based on Hankel factors and horizontal
tx-rx distance. Shape: (Hankel, n_rx)
model: np.ndarray
Polar coords of the receiver pos (3, n).
rho: np.ndarray
Resistivities for each layer (Ohm*m).
d: np.ndarray
Thicknesses of each layer (m).
f: float
Frequency (Hz).
mode: str
Calculation for 'te', 'tm' or 'tetm' possible. Mode 'te' for closed
loops and 'tetm' for grounded wires. Single 'tm' is only for debug.
Returns
-------
aa : np.ndarray
Ratio of the partial wave amplitude A(z,u)/A(0,u)
aap : np.andarray
Ratio of the partial wave amplitude A'(z,u)/A'(0,u)
bt : np.ndarray
Admittance at the surface of the layerd halfspace
"""
aa = np.ones(u.shape, dtype=np.complex) #
aap = np.ones(u.shape, dtype=np.complex) # shape = hankel x model
bt = np.ones(u.shape, dtype=np.complex) #
z_idx_air = np.where(model[2] >= 0.)[0]
z_idx_halfspace = np.where(model[2] < 0.)[0]
if len(z_idx_air) != 0:
bt[:, z_idx_air] = btp(u[:, z_idx_air], model[:, z_idx_air],
rho, d, f, mode)
if len(z_idx_halfspace) != 0:
bt[:, z_idx_halfspace], aa[:, z_idx_halfspace], \
aap[:, z_idx_halfspace] =\
downward(u[:, z_idx_halfspace],
model[:, z_idx_halfspace], rho, d, f, mode)
return aa, aap, bt
[docs]def hfield_3D_hed_te(polar, u, aa, aap, bt, f, Ids):
""" Calculation of magnetic field for transversal electric mode.
Computes the transversal electric induced magnetic field of a x-directed
dipole at (0, 0, 0). Field shape (3, n) with x, y, z components for each
reciever point in *polar*.
Internal function. Called by **makeField** if *ftype* == 'H' and *mode* in
('te', 'tetm').
Parameters
----------
polar: np.ndarray
Polar coords of the receiver pos (3, n).
u: np.ndarray
Horizontal wavenumbers based on Hankel factors and horizontal
tx-rx distance. Shape: (Hankel, n_rx)
aa : np.ndarray
Ratio of the partial wave amplitude A(z,u)/A(0,u)
aap : np.andarray
Ratio of the partial wave amplitude A'(z,u)/A'(0,u)
bt : np.ndarray
Admittance at the surface of the layerd halfspace
f: float
Frequency (Hz).
Ids: float
Dipole current * dipole length.
Returns
-------
field: np.ndarray
Transversal electric induced magnetic field of a x-directed dipole
at (0, 0, 0). field.shape = polar.shape.
"""
# empty field vector with dim=dim(polar)
field = np.zeros(polar.shape, dtype=np.complex)
# hankel filter coefficients
fc0, nc, nc0 = hankelfc(3)
fc1, nc, nc0 = hankelfc(4)
# indices of the meshpoints in the airspace
idx_air = np.array(np.nonzero(polar[2] >= 0.))
if len(idx_air.shape) > 1:
idx_air = idx_air[0]
# indices of the meshpoints inside the layered earth
idx_layer = np.array(np.nonzero(polar[2] < 0.))
if len(idx_layer.shape) > 1:
idx_layer = idx_layer[0]
# core function
delta = 2./(u + bt)
delta[:, idx_air] *= np.exp(u[:, idx_air] * polar[2][idx_air] * (-1.))
# convolution
aux1 = np.zeros((2, polar.shape[1]), dtype=np.complex)
aux2 = np.zeros((1, polar.shape[1]), dtype=np.complex)
# the six hankel integrals ! 3 for air, 3 for the layered earth
aux1[0, idx_air] = np.sum((delta[:, idx_air] *
u[:, idx_air] *
fc1[::-1]),
0, np.complex)
aux1[0, idx_layer] = np.sum((delta[:, idx_layer] *
bt[:, idx_layer] *
aap[:, idx_layer] *
fc1[::-1]),
0, np.complex)
aux1[1, idx_air] = np.sum((delta[:, idx_air] *
u[:, idx_air] *
u[:, idx_air] *
fc1[::-1]),
0, np.complex)
aux1[1, idx_layer] = np.sum((delta[:, idx_layer] *
u[:, idx_layer] *
u[:, idx_layer] *
aa[:, idx_layer] *
fc1[::-1]),
0, np.complex)
aux2[0, idx_air] = np.sum((delta[:, idx_air] *
u[:, idx_air] *
u[:, idx_air] *
fc0[::-1]),
0, np.complex)
aux2[0, idx_layer] = np.sum((delta[:, idx_layer] *
u[:, idx_layer] *
bt[:, idx_layer] *
aap[:, idx_layer] *
fc0[::-1]),
0, np.complex)
aux1 /= polar[0]
aux2 /= polar[0]
# field
bh = 1./(4.*np.pi)
sphi = np.sin(polar[1])
cphi = np.cos(polar[1])
field[2] = Ids * bh * sphi * aux1[1]
# only for r, phi components
sphi[idx_layer] *= -1.
cphi[idx_layer] *= -1.
field[0] = -Ids * bh * sphi * 1./polar[0] * aux1[0] + \
Ids * bh * sphi * aux2[0]
field[1] = Ids * bh * cphi * 1./polar[0] * aux1[0]
field_cartesian = vec.PtoK_field(polar, field, dtype=np.complex128)
# change the direction of the field due to the z-convention
field_cartesian[0] *= -1
field_cartesian[1] *= -1
# field_cartesian[0] = np.conjugate(field_cartesian[0]) # 22.05 -x_imag
return field_cartesian
[docs]def hfield_3D_hed_tm(polar, u, aa, aap, bt, f, Ids):
""" Calculation of magnetic field for transversal magnetic mode.
Computes the transversal magnetic component of the magnetic field of a
x-directed dipole at (0, 0, 0). Field shape (3, n) with x, y, z components
for each reciever point in *polar*.
Internal function. Called by **makeField** if *ftype* == 'H' and *mode* in
('tm', 'tetm').
Parameters
----------
polar: np.ndarray
Polar coords of the receiver pos (3, n).
u: np.ndarray
Horizontal wavenumbers based on Hankel factors and horizontal
tx-rx distance. Shape: (Hankel, n_rx)
aa : np.ndarray
Ratio of the partial wave amplitude A(z,u)/A(0,u)
aap : np.andarray
Ratio of the partial wave amplitude A'(z,u)/A'(0,u)
bt : np.ndarray
Admittance at the surface of the layerd halfspace
f: float
Frequency (Hz).
Ids: float
Dipole current * dipole length.
Returns
-------
field: np.ndarray
Transversal magnetic component of the magnetic field of a
x-directed dipole at (0, 0, 0). field.shape = polar.shape.
"""
# empty field vector with dim=dim(polar)
field = np.zeros(polar.shape, dtype=np.complex)
# hankel filter coefficients
fc0, nc, nc0 = hankelfc(3)
fc1, nc, nc0 = hankelfc(4)
# indices of the meshpoints in the airspace
idx_air = np.array(np.nonzero(polar[2] >= 0.))
if len(idx_air.shape) > 1:
idx_air = idx_air[0]
# indices of the meshpoints inside the layered earth
idx_layer = np.array(np.nonzero(polar[2] < 0.))
if len(idx_layer.shape) > 1:
idx_layer = idx_layer[0]
# core function
delta = np.zeros(np.shape(u))
# tm air == 0 !
delta[:, idx_layer] = 2.0
# convolution
aux1 = np.zeros((polar.shape[1]), dtype=np.complex)
aux2 = np.zeros((polar.shape[1]), dtype=np.complex)
# the two hankel integrals ! 2 for air, 2 for the layered earth
aux1[idx_air] = np.sum((delta[:, idx_air] *
fc1[::-1]),
0, np.complex)
aux1[idx_layer] = np.sum((delta[:, idx_layer] *
aa[:, idx_layer] * # aa is valid
fc1[::-1]),
0, np.complex)
aux2[idx_air] = np.sum((delta[:, idx_air] *
u[:, idx_air] *
fc0[::-1]),
0, np.complex)
aux2[idx_layer] = np.sum((delta[:, idx_layer] *
u[:, idx_layer] *
aa[:, idx_layer] * # aa is valid
fc0[::-1]),
0, np.complex)
aux1 /= polar[0]
aux2 /= polar[0]
# field
bh = 1./(4.*np.pi)
sphi = np.sin(polar[1])
cphi = np.cos(polar[1])
field[0] = -Ids * bh * sphi * (1./polar[0]) * aux1
field[1] = (Ids * bh * cphi * (1./polar[0]) * aux1) - \
(Ids * bh * cphi * aux2)
field_cartesian = vec.PtoK_field(polar, field, dtype=np.complex128)
# change the direction of the field due to the z-convention
field_cartesian[0] *= -1
field_cartesian[1] *= -1
return field_cartesian
[docs]def efield_3D_hed_te(polar, u, aa, aap, bt, f, Ids):
""" Calculation of electric field for transversal electric mode.
Computes the transversal electric induced electric field of a
x-directed dipole at (0, 0, 0). Field shape (3, n) with x, y, z components
for each reciever point in *polar*.
Internal function. Called by **makeField** if *ftype* == 'E' and *mode* in
('te', 'tetm').
Parameters
----------
polar: np.ndarray
Polar coords of the receiver pos (3, n).
u: np.ndarray
Horizontal wavenumbers based on Hankel factors and horizontal
tx-rx distance. Shape: (Hankel, n_rx)
aa : np.ndarray
Ratio of the partial wave amplitude A(z,u)/A(0,u)
aap : np.andarray
Ratio of the partial wave amplitude A'(z,u)/A'(0,u)
bt : np.ndarray
Admittance at the surface of the layerd halfspace
f: float
Frequency (Hz).
Ids: float
Dipole current * dipole length.
Returns
-------
field: np.ndarray
Transversal electric component of the electric field of a
x-directed dipole at (0, 0, 0). field.shape = polar.shape.
"""
# empty field vector with dim=dim(polar)
field = np.zeros(polar.shape, dtype=np.complex)
# hankel filter coefficients
fc0, nc, nc0 = hankelfc(3)
fc1, nc, nc0 = hankelfc(4)
# indices of the meshpoints in the airspace
idx_air = np.array(np.nonzero(polar[2] >= 0.))
if len(idx_air.shape) > 1:
idx_air = idx_air[0]
# indices of the meshpoints inside the layered earth
idx_layer = np.array(np.nonzero(polar[2] < 0.))
if len(idx_layer.shape) > 1:
idx_layer = idx_layer[0]
# core function:
# air: 2 / (lam + u1) * exp(lam*z)
# ground: 2 / (lam + u1)
delta = 2./(u + bt)
delta[:, idx_air] *= np.exp(u[:, idx_air] * polar[2][idx_air] * (-1.))
# convolution
aux1 = np.zeros((polar.shape[1]), dtype=np.complex)
aux2 = np.zeros((polar.shape[1]), dtype=np.complex)
# 2 / (lam + u1) * exp(lam*z) * J1
aux1[idx_air] = np.sum((delta[:, idx_air] *
fc1[::-1]),
0, np.complex)
# 2 / (lam + u1) * exp(-u1*z) * J1
aux1[idx_layer] = np.sum((delta[:, idx_layer] *
aa[:, idx_layer] *
fc1[::-1]),
0, np.complex)
aux2[idx_air] = np.sum((delta[:, idx_air] *
u[:, idx_air] *
fc0[::-1]),
0, np.complex)
aux2[idx_layer] = np.sum((delta[:, idx_layer] *
u[:, idx_layer] *
aa[:, idx_layer] *
fc0[::-1]),
0, np.complex)
aux1 /= polar[0] # bessel/hankel
aux2 /= polar[0] # bessel/hankel
# field
bh = 1./(4.*np.pi)
sphi = np.sin(polar[1])
cphi = np.cos(polar[1])
mu0 = 4e-7 * np.pi # mu null
iwm = 1j * 2. * np.pi * f * mu0
field[0] = -Ids * iwm * bh * cphi * 1./polar[0] * aux1
field[1] = -Ids * iwm * bh * sphi * 1./polar[0] * aux1 + \
Ids * iwm * bh * sphi * aux2
field_cartesian = vec.PtoK_field(polar, field, dtype=np.complex128)
# change the direction of the field due to the z-convention
# e^-iwt conformity with custEM: added 19.11.2019 *= -1 to all
# field_cartesian[0] *= -1
# field_cartesian[1] *= -1
field_cartesian[2] *= -1
return field_cartesian
[docs]def calcField(polar, rho, d, f, Ids, ftype, mode):
""" Calculates field for a given dipole on given polar coords. Internally.
Internally used by *makeField*. Please be referred to the docstrings of
**makeField**. And please use makeField directly!!!
"""
# calculate horizontal wavenumbers
fc0, nc, nc0 = hankelfc(3)
nu = np.arange(0, nc, 1, dtype=np.int) + 1
n = nc0 - nc + nu
# n = [-39, ...., 60]
q = 0.1 * np.log(10)
# q = 0.23025851
# 23.02... samples per decade (Anderson (1982), pp.364)
hlp = np.exp(-(n - 1) * q)
u = hlp[:, np.newaxis] / polar[0] # one u per receiver coord
RHO = np.copy(rho)
D = np.copy(d)
if mode == 'te': # closed loops
aaTE, aapTE, btTE = downout(u, polar, RHO, D, f, 'te')
if ftype.upper() == 'E':
F = efield_3D_hed_te(polar, u, aaTE, aapTE, btTE, f, Ids)
elif ftype.upper() in ('B', 'H'):
F = hfield_3D_hed_te(polar, u, aaTE, aapTE, btTE, f, Ids)
elif mode == 'tm': # testing purpose
aaTM, aapTM, btTM = downout(u, polar, RHO, D, f, 'tm')
if ftype.upper() == 'E':
raise Exception('E + tm: This should not happen. This '
'should be prevented by the function makeField.')
elif ftype.upper() in ('B', 'H'):
F = hfield_3D_hed_tm(polar, u, aaTM, aapTM, btTM, f, Ids)
elif mode == 'tetm': # grounded wires
aaTE, aapTE, btTE = downout(u, polar, RHO, D, f, 'te')
aaTM, aapTM, btTM = downout(u, polar, RHO, D, f, 'tm')
if ftype.upper() == 'E':
raise Exception('E + tm: This should not happen. This '
'should be prevented by the function makeField.')
elif ftype.upper() in ('B', 'H'):
F = hfield_3D_hed_tm(polar, u, aaTM, aapTM, btTM, f, Ids)
F += hfield_3D_hed_te(polar, u, aaTE, aapTE, btTE, f, Ids)
if ftype.upper() == 'B':
F *= 4e-7 * np.pi # mu null
return F
[docs]def makeField(coords, rho_in, d_in, f=2000, Ids=1, pos=(0, 0, 0), angle=0,
mode='te', inputType='M', ftype='B',
cell_center=False, drop_tol=1e-2, src_z=-0.01):
'''
Calculation of the electric or magnetic field of a horizontal electric
dipole at position pos, pointing in a direction defined by angle on
given cartesian coordinates.
Parameters
----------
coords: np.ndarray or string
Reciever coords. Possible input types are numpy ndarrays for direct
cartesian coordinates, ranges for (x, y, z) or pygimli Meshes.
rho_in: float or np.ndarray
Float or Array of resistivity values for the 1d layered earth model.
Airspace is at the level of the source dipole.
d_in: float or np.ndarray
Layer thicknesses in m. As the lower halfspace is considered to have an
infinite thickness, *d_in* is always one value short of *rho_in* (an
empty list ar array or a 0 for homogenous halfspace.)
f: float [ 2000 ]
Frequency (Hz).
Ids: float [ 1 ]
Dipole current * dipole length. Used for simple scaling of the
calcualted field.
pos: tuple of length 3 [ (0 ,0, 0) ]
Absolute position of the source dipole in cartesian coordinates.
Values for z are used for a shift of the airspace. Currently only
sources at the upper halfpsce boundary are permitted.
angle: float [ 0 ]
Rotation of the dipole with respect of an x-directed dipole counting
positive clockwise.
mode: string [ 'te' ]
For a closed loop consisting of a finite number of dipoles the total
field can be seen as superposition of the transversal electric
components of the single dipole fields ('te'). For grounded dipoles
'tetm' is needed.
inputType: string [ 'M' ]
Specifier for input coordinates. Possible choices are 'M' if *coords*
is a pygimli mesh or file path to a pygimli mesh, 'C' if coords is
a np.ndarray with ranges to build a meshgrid, or 'V' to indicate that
*coords* is a vector of cartesian coordinates.
ftype: string [ 'B' ]
Flag to control calculated field type. Possible choices are 'E', 'B' or
'H' (asuming B = 4e-7 * pi * H).
verbose: boolean [ False ]
Turn on verbose mode.
cell_center: boolean [ False ]
If coords is a pygimli mesh, there is the additional possibility to
calulate the fields in the cell Centers, instead of the node positions.
drop_tol: float [ 1e-2 ]
Singularity fix. All horizontal distances between *drop_tol* and the
transmitter dipole are placed between the first reciever outside the
tolerance and the tolerance, maintaining the correct order and angle.
This has been very useful for later usage of the fields in FEM
approaches.
src_z: float [-0.001]
This is only used if grounded terms for an electric field are used. In
this case the source has to be buried in order to get the correct
results. Default is 1 cm (remember: z defined positive upwards). So in
most cases this value should be negative.
'''
# %% handle receiver coords
transposed = False
if inputType == 'C':
x = coords[0]
y = coords[1]
z = coords[2]
grid = np.meshgrid(x, y, z, indexing='ij')
modelvector = vec.GridtoVector(grid)
elif inputType == 'V':
if np.shape(coords)[1] == 3 and np.shape(coords)[0] != 3:
modelvector = coords.T
transposed = True
else:
modelvector = coords
elif inputType == 'M':
# no isinstance(coords, pg.Mesh) due to independance of pygimli
if not cell_center:
modelvector = coords.positions().array().T
else:
modelvector = coords.cellCenters().array().T
elif inputType == 'Grid':
# only used in some test scripts, therefore not mentioned in the
# docstring
modelvector = vec.GridtoVector(coords)
# not origin == (0, 0, 0)
if not np.allclose(pos, np.array((0, 0, 0))):
model = vec.translate(modelvector, -pos[0], -pos[1], z=-pos[2],
copy=True)
else:
model = np.copy(modelvector)
# handle resistivity distribution
rho = np.atleast_1d(rho_in)
d = np.atleast_1d(d_in)
if len(d) == 1 and np.isclose(d[0], 0):
d = np.array([])
# %% field calculation started
if ftype.upper() == 'E' and mode != 'te':
# e field tetm
e = World1D(rho=rho, thk=d, f=f)
hed = HED(src_z=src_z, src_theta=angle, src_ids=Ids,
config=e, timer=None, debug=False)
hed.setCoords(model)
hed.calculate()
F = hed.e_field
else:
# cartesian -> cylindrical coordinates
cylindrical = vec.KtoP(model, drop_tol=drop_tol)
if not np.isclose(angle, 0):
cylindrical[1] -= angle # rotation of the polar coords
F = calcField(cylindrical, rho, d, f, Ids, ftype, mode)
if not np.isclose(angle, 0):
vec.rotate3(F, angle) # rotation of field in cartesian coords
# %% handling phi (rotate field back) + return
# return field in correct shape
if inputType == 'Grid':
return vec.VectortoGrid(F, np.shape(coords)[1:])
if transposed:
return F.T
else:
return F
# The End