"""
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/>.
from comet.pyhed.misc import fixSingularity
import numpy as np
[docs]def hedx_magnetic(model, f, sigma, I, ds, drop_tol=1e-6):
"""
Analytic calculation of the magnatic field for an electric dipole in x
direction. Formula given in Ward and Hohmann (1988), page 174 number 2.41.
Model in cartesian coordinates, as well as the output.
Sigma != 0
drop_tol to avoid singularities [1e-6]
No grounding!
"""
assert not np.isclose(sigma, 0.0), 'The function won\'t work for vakuum \
space. The formulations are not suited for the resistive case.'
radius = np.sqrt(model[0]**2 + model[1]**2 + model[2]**2)
model = fixSingularity(model, drop_tol=drop_tol, dipole_z=0.0)
k = np.sqrt((0 + 1j) * 2 * np.pi * f * 4e-7 * np.pi * sigma)
ikr = 1j * k * radius
factor = ((I * ds) / (4 * np.pi * radius**2)) *\
(ikr + 1) * np.exp(-ikr)
hy = factor * (-model[2] / radius)
hz = factor * (model[1] / radius)
hx = np.zeros_like(hy, complex)
return np.array((hx, hy, hz))
[docs]def hedx_electric(model, f, sigma, I, ds, drop_tol=1e-6):
"""
Analytic calculation of the electric field for an electric dipole in x
direction. Formula given in Ward and Hohmann (1988), page 173 number 2.40.
Model in cartesian coordinates, as well as the output.
Sigma != 0
drop_tol to avoid singularities [1e-6]
No grounding!
"""
assert not np.isclose(sigma, 0), 'The function won\'t work for vakuum \
space. The formulations are not suited for the resistive case.'
radius = np.sqrt(model[0]**2 + model[1]**2 + model[2]**2)
radius = fixSingularity(radius, drop_tol=drop_tol)
k = np.sqrt((0 + 1j) * 2 * np.pi * f * 4e-7 * np.pi * sigma)
ikr = 1j * k * radius
kkrr = (radius**2) * (k**2)
factor = ((I * ds) / (4 * np.pi * sigma * radius**3)) * np.exp(-ikr)
ex = factor * (((model[0] / radius)**2) * (-kkrr + 3 * ikr + 3) +
(kkrr - ikr - 1))
ey = factor * (model[0] * model[1] / (radius**2)) * (-kkrr + 3 * ikr + 3)
ez = factor * (model[0] * model[2] / (radius**2)) * (-kkrr + 3 * ikr + 3)
return np.array((ex, ey, ez))
# The End