"""
Part of comet/pyhed/hed/reference
"""
# 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/>.
import numpy as np
import scipy.special as scsp
import sys
from comet.pyhed.misc import KtoP
[docs]def hed_field(r, f, sigma, phi, I, ds, BorH='B'):
''' Semi analytic solution for electrical and magnetical fields at the
surface of a homogeneous halfspace of conductivity sigma.
'''
phi *= -1 # temporary change
k = np.sqrt((0 + 1j) * 2 * np.pi * f * 4e-7 * np.pi * sigma)
u = k * r
u2 = u / 2
# phi += np.pi / 2
# bessel factors
I0 = scsp.iv(0, u2) # first kind, order 0
I1 = scsp.iv(1, u2) # first kind, order 1
K0 = scsp.kv(0, u2) # second kind, order 0
K1 = scsp.kv(1, u2) # second kind, order 1
# exponential term
E = np.exp(-u)
# coefficients for the magnetic field components
co1 = I * ds * np.sin(phi) / (2 * np.pi * (r**2))
co2 = I * ds * np.cos(phi) / (2 * np.pi * (r**2))
# coefficients for the electric field components
co3 = I * ds * np.sin(phi) / (2 * np.pi * sigma * (r ** 3))
co4 = I * ds * np.cos(phi) / (2 * np.pi * sigma * (r ** 3))
# electric field
Er = 1 + (1 + u) * E
Ef = 2 - (1 + u) * E
Ez = -(u ** 2) * I1 * K1
er = co4 * Er
ef = co3 * Ef
ez = co4 * Ez
# magnetic field
if BorH == 'B':
borh = 1. # A/m
elif BorH == 'H':
borh = 4 * np.pi * 100 # nT
Hr = -(u2) * ((I1 * K0 - I0 * K1) - 3 * I1 * K1)
Hf = I1 * K1
Hz = (3 - (3 + 3 * u + (u ** 2)) * E) / (u ** 2)
# -(u**2) im ersten Term laut Hohmann
# Hz = (3 - (3 + 3 * u - (u ** 2)) * E) / (u ** 2)
hr = borh * co1 * Hr
hf = borh * co2 * Hf
hz = borh * co1 * Hz
return (hr, hf, hz), (er, ef, ez) # , (I0, I1, K0, K1, u)
[docs]def hed_field_hohmann(model, f, sigma, I, ds, ftype='H', **kwargs):
'''
semi analytic solution for a magnetic field at the
surface of a homogeneous halfspace of conductivity sigma
ward hohmann formula: page 235-236 No 4.166, 4.171 and 4.173
edit:
E-term: grounding term only (4.159)
'''
# radial distance
# model[1] *= -1
polar = KtoP(model, drop_tol=kwargs.pop('drop_tol', 1e-6))
k = np.sqrt(-1j * 2 * np.pi * f * 4e-7 * np.pi * sigma)
u = 1j * k * polar[0] # ikp
u2 = u / 2
# bessel factors
I0 = scsp.iv(0, u2) # first kind, order 0
I1 = scsp.iv(1, u2) # first kind, order 1
K0 = scsp.kv(0, u2) # second kind, order 0
K1 = scsp.kv(1, u2) # second kind, order 1
# exponential term
E = np.exp(-u)
if ftype == 'E':
# grounding term !
Ex = (I * ds / (2 * np.pi * sigma * (polar[0]**3))) * \
(-2 + (u + 1)*E - (3 * (model[0]**2) / (polar[0]**2)))
Ey = (I * ds / (2 * np.pi * sigma * (polar[0]**3))) * \
(1 + (u + 1)*E - (3 * (model[1]**2) / (polar[0]**2)))
# Ex = -(I * ds / (2 * np.pi * sigma * (polar[0]**3))) * \
# (-1 + (u + 1)*E + (model[1]**2 - 2 * model[0]**2) / (polar[0]**2))
# Ey = (I * ds / (2 * np.pi * sigma * (polar[0]**6))) * \
# (3 * model[0] * model[1])
return np.array((Ex, Ey, np.zeros_like(Ey)),
np.complex128)
Hx = ((I * ds * model[0] * model[1]) / (4 * np.pi * (polar[0]**4))) * \
(u * (I0 * K1 - I1 * K0) - 8 * I1 * K1)
Hy = ((-I * ds) / (4 * np.pi * (polar[0]**2))) * \
(6 * I1 * K1 + u * (I1 * K0 - I0 * K1) + ((model[0]**2) /
(polar[0]**2)) * (u * (I0 * K1 - I1 * K0) - 8 * I1 * K1))
Hz = ((-I * ds * model[1]) / (2 * np.pi * (k**2) * (polar[0]**5))) * \
(3 - (3 + 3 * u + (u**2)) * E)
return np.array((Hx, Hy, Hz), np.complex128)
if __name__ == '__main__':
hed_field(sys.argv[1:])
# THE END