# Author: LKouadio <etanoyau@gmail.com>
# License: LGPL-3.0
"""
TEM time-domain → frequency-domain transforms.
Three loop configurations are fully supported:
**Central / coincident loop** (``offset = 0``)
Receiver at the centre of the transmitter loop. Uses the Ward &
Hohmann (1988) late-time formula directly:
.. math::
\\rho_a(t) = \\left(
\\frac{M\\,\\mu_0^{5/2}}
{10\\sqrt{\\pi}\\,|\\partial B_z/\\partial t|\\,t^{5/2}}
\\right)^{2/3}
**In-loop / large-loop** (``0 < offset < inner_radius``)
Receiver is inside the transmitter loop but not at the centre.
A Biot-Savart static-field correction factor
:math:`\\eta = H_z^{\\rm static}(r_x, r_y)\\,/\\,H_z^{\\rm static}(0,0)`
is applied to the effective transmitter moment:
.. math::
\\rho_a(t) = \\left(
\\frac{M\\,\\eta\\,\\mu_0^{5/2}}
{10\\sqrt{\\pi}\\,|\\partial B_z/\\partial t|\\,t^{5/2}}
\\right)^{2/3}
The Biot-Savart integral is evaluated analytically for rectangular /
square loops and numerically for circular loops.
**Offset / separated loop** (``offset ≥ inner_radius``)
Transmitter and receiver are separated by horizontal distance *d*.
Uses the Ward & Hohmann (1988) offset-dipole formula:
.. math::
\\rho_a(t) = \\left(
\\frac{M\\,\\mu_0^{5/2}}
{20\\sqrt{\\pi}\\,|\\partial B_z/\\partial t|\\,d^3\\,t^{5/2}}
\\right)^{2/3}
:class:`LateTimeTransform`
Fast, approximate method — standard industry approach.
:class:`FourierTransform`
Rigorous Fourier cosine transform with Kramers-Kronig reconstruction.
Valid at all time gates; first-order waveform deconvolution for all
waveform types (Fitterman & Stewart 1986).
:class:`TEMtoEDI`
High-level dispatcher → :class:`~pycsamt.seg.collection.EDICollection`.
"""
from __future__ import annotations
from collections.abc import Sequence
import numpy as np
from ..api.property import PyCSAMTObject
from ._base import TEMSounding
__all__ = [
"LateTimeTransform",
"FourierTransform",
"TEMtoEDI",
# geometry helpers (importable for testing)
"build_adjacency",
"_rho_a_late_time",
"_rho_a_in_loop",
"_rho_a_offset_loop",
"_biot_savart_rect_hz",
"_biot_savart_circle_hz",
"_biot_savart_rect_hz_segments",
"_biot_savart_circle_hz_segments",
"_in_loop_geometry_factor",
"_in_loop_geometry_factor_td",
"_pseudo_freq",
"_phase_from_rho",
"_build_z_array",
"_z_error_from_rho_error",
"_cosine_transform_1d",
"_kramers_kronig_re",
"_waveform_moments",
"_apply_waveform_correction",
"MU0",
]
MU0 = 4.0 * np.pi * 1e-7 # H/m
# ---------------------------------------------------------------------------
# Biot-Savart static field helpers
# ---------------------------------------------------------------------------
def _hz_horiz_seg(
y0: float, x1: float, x2: float, rx: float, ry: float
) -> float:
"""Hz from a horizontal wire segment at y=y0, running x1→x2."""
d = ry - y0
if abs(d) < 1e-14:
return 0.0
u1, u2 = x1 - rx, x2 - rx
return (1.0 / d) * (
u2 / np.sqrt(d * d + u2 * u2) - u1 / np.sqrt(d * d + u1 * u1)
)
def _hz_vert_seg(
x0: float, y1: float, y2: float, rx: float, ry: float
) -> float:
"""Hz from a vertical wire segment at x=x0, running y1→y2."""
d = rx - x0
if abs(d) < 1e-14:
return 0.0
u1, u2 = y1 - ry, y2 - ry
return -(1.0 / d) * (
u2 / np.sqrt(d * d + u2 * u2) - u1 / np.sqrt(d * d + u1 * u1)
)
def _biot_savart_rect_hz(
rx: float,
ry: float,
a: float,
b: float,
) -> float:
r"""
Static :math:`H_z` per unit current at ``(rx, ry)`` inside a
rectangular transmitter loop of half-sides ``(a, b)`` centred at
the origin (analytical Biot-Savart).
The loop current runs counter-clockwise when viewed from above:
bottom (+x) → right (+y) → top (−x) → left (−y).
Parameters
----------
rx, ry : float
Receiver position in metres. Must be strictly inside
the loop (|rx| < a, |ry| < b) to avoid a singularity.
a, b : float
Half-sides of the rectangular loop (m). For a square loop
pass ``a = b = side / 2``.
Returns
-------
float
:math:`H_z` in A m⁻¹ per ampere of transmitter current.
References
----------
.. [1] Griffiths, D.J. (1999). *Introduction to Electrodynamics*,
3rd ed., Prentice Hall — §5.4 (Biot-Savart law for line current).
"""
hz = (
_hz_horiz_seg(-b, -a, +a, rx, ry) # bottom: (−a,−b) → (+a,−b)
+ _hz_vert_seg(+a, -b, +b, rx, ry) # right: (+a,−b) → (+a,+b)
+ _hz_horiz_seg(+b, +a, -a, rx, ry) # top: (+a,+b) → (−a,+b)
+ _hz_vert_seg(-a, +b, -b, rx, ry) # left: (−a,+b) → (−a,−b)
)
return hz / (4.0 * np.pi)
def _biot_savart_rect_hz_segments(
rx: float,
ry: float,
a: float,
b: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
r"""
Per-segment :math:`H_z / (4\pi)` and perpendicular distances for a
rectangular transmitter loop of half-sides ``(a, b)``.
Returns four arrays of length 4 (segments: bottom, right, top, left):
* ``hz_rx`` — Hz contribution at ``(rx, ry)`` per segment.
* ``hz_00`` — Hz contribution at ``(0, 0)`` per segment.
* ``dist_rx`` — Perpendicular distance from ``(rx, ry)`` to each
segment's line.
* ``dist_00`` — Perpendicular distance from ``(0, 0)`` to each
segment's line.
Used by :func:`_in_loop_geometry_factor_td` to weight each segment's
contribution by the EM diffusion length at each time gate.
"""
hz_rx = np.array(
[
_hz_horiz_seg(-b, -a, +a, rx, ry), # bottom
_hz_vert_seg(+a, -b, +b, rx, ry), # right
_hz_horiz_seg(+b, +a, -a, rx, ry), # top
_hz_vert_seg(-a, +b, -b, rx, ry), # left
]
) / (4.0 * np.pi)
hz_00 = np.array(
[
_hz_horiz_seg(-b, -a, +a, 0.0, 0.0),
_hz_vert_seg(+a, -b, +b, 0.0, 0.0),
_hz_horiz_seg(+b, +a, -a, 0.0, 0.0),
_hz_vert_seg(-a, +b, -b, 0.0, 0.0),
]
) / (4.0 * np.pi)
# Perpendicular distances from (rx,ry) and (0,0) to each segment's line
dist_rx = np.array([abs(ry + b), abs(a - rx), abs(b - ry), abs(rx + a)])
dist_00 = np.array([b, a, b, a])
return hz_rx, hz_00, dist_rx, dist_00
def _biot_savart_circle_hz(
rx: float,
ry: float,
radius: float,
n_seg: int = 720,
) -> float:
r"""
Static :math:`H_z` per unit current at ``(rx, ry)`` inside a
circular transmitter loop of given ``radius`` (numerical
Biot-Savart via piecewise-linear integration).
Parameters
----------
rx, ry : float
Receiver position (m). Should be inside the loop (r < radius).
radius : float
Loop radius (m).
n_seg : int
Number of equal arc segments for the integration. 720 gives
relative errors < 10⁻⁵ for rx/radius < 0.9.
Returns
-------
float
:math:`H_z` in A m⁻¹ per ampere of transmitter current.
"""
theta = np.linspace(0.0, 2.0 * np.pi, n_seg, endpoint=False)
dtheta = 2.0 * np.pi / n_seg
xw = radius * np.cos(theta)
yw = radius * np.sin(theta)
dlx = -radius * np.sin(theta) * dtheta # tangent × arc length
dly = +radius * np.cos(theta) * dtheta
dx = rx - xw
dy = ry - yw
r3 = (dx * dx + dy * dy) ** 1.5
hz = np.sum((dlx * dy - dly * dx) / r3)
return hz / (4.0 * np.pi)
def _biot_savart_circle_hz_segments(
rx: float,
ry: float,
radius: float,
n_seg: int = 360,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
r"""
Per-element :math:`H_z / (4\pi)` and distances for a circular loop.
The circle is discretised into ``n_seg`` equal arc elements. For each
element the Biot-Savart contribution at ``(rx, ry)`` and at ``(0, 0)``
is returned together with the element-to-receiver distances.
Returns
-------
hz_rx : ndarray (n_seg,)
hz_00 : ndarray (n_seg,)
dist_rx : ndarray (n_seg,)
Distance from each element to ``(rx, ry)``.
dist_00 : ndarray (n_seg,)
Distance from each element to ``(0, 0)`` = ``radius`` (uniform).
"""
theta = np.linspace(0.0, 2.0 * np.pi, n_seg, endpoint=False)
dt = 2.0 * np.pi / n_seg
xw = radius * np.cos(theta)
yw = radius * np.sin(theta)
dlx = -radius * np.sin(theta) * dt
dly = +radius * np.cos(theta) * dt
# At (rx, ry)
dx_rx = rx - xw
dy_rx = ry - yw
r3_rx = (dx_rx**2 + dy_rx**2) ** 1.5
hz_rx = (dlx * dy_rx - dly * dx_rx) / r3_rx / (4.0 * np.pi)
# At (0, 0) — all wire elements at the same distance = radius
hz_00 = (dlx * (-yw) - dly * (-xw)) / radius**3 / (4.0 * np.pi)
dist_rx = np.sqrt(dx_rx**2 + dy_rx**2)
dist_00 = radius * np.ones(n_seg)
return hz_rx, hz_00, dist_rx, dist_00
def _in_loop_geometry_factor(
rx: float,
ry: float,
loop_shape: str,
loop_dims: tuple,
) -> float:
r"""
Geometric correction factor
:math:`\eta = H_z^{\rm static}(r_x, r_y)\,/\,H_z^{\rm static}(0, 0)`
for a receiver inside the transmitter loop at position ``(rx, ry)``
relative to the loop centre.
Used by :func:`_rho_a_in_loop` to correct the effective transmitter
moment for non-central in-loop receiver positions.
Parameters
----------
rx, ry : float
Receiver offset from loop centre (m).
loop_shape : str
``'square'``, ``'rectangular'``, or ``'circle'``.
loop_dims : tuple of float
Shape parameters:
* ``'square'`` → ``(side_length,)``
* ``'rectangular'`` → ``(length, width)``
* ``'circle'`` → ``(radius,)``
Returns
-------
float
Correction factor η > 0. Equal to 1.0 when (rx, ry) = (0, 0).
Notes
-----
At late time (diffusion time >> loop size), η has no effect because
the induced eddy current "smoke-ring" has expanded far beyond the
loop and the field is spatially uniform. The correction is most
significant at early-to-mid time gates and for large offsets r/L.
References
----------
.. [1] Nabighian, M.N. & Macnae, J.C. (1991). Time domain
electromagnetic prospecting methods. In *Electromagnetic Methods
in Applied Geophysics*, Vol. 2, pp. 427–520. SEG.
.. [2] Spies, B.R. & Frischknecht, F.C. (1991). Electromagnetic
sounding. Ibid., pp. 285–425.
"""
if abs(rx) < 1e-10 and abs(ry) < 1e-10:
return 1.0
if loop_shape in ("square", "rectangular"):
if loop_shape == "square":
a = b = float(loop_dims[0]) / 2.0
else:
a = float(loop_dims[0]) / 2.0
b = float(loop_dims[1]) / 2.0
hz_rx = _biot_savart_rect_hz(rx, ry, a, b)
hz_00 = _biot_savart_rect_hz(0.0, 0.0, a, b)
elif loop_shape == "circle":
r = float(loop_dims[0])
hz_rx = _biot_savart_circle_hz(rx, ry, r)
hz_00 = _biot_savart_circle_hz(0.0, 0.0, r)
else:
# Unknown shape — treat as equivalent-area square
a = b = np.sqrt(float(loop_dims[0])) / 2.0
hz_rx = _biot_savart_rect_hz(rx, ry, a, b)
hz_00 = _biot_savart_rect_hz(0.0, 0.0, a, b)
if hz_00 < 1e-30:
return 1.0
return float(max(hz_rx / hz_00, 0.01)) # guard against unphysical values
def _in_loop_geometry_factor_td(
rx: float,
ry: float,
loop_shape: str,
loop_dims: tuple,
t: np.ndarray,
rho_a: np.ndarray,
) -> np.ndarray:
r"""
Time-dependent geometry factor :math:`\eta(t)` for in-loop receivers.
Replaces the DC Biot-Savart approximation with a physically correct
transition between the near-field (early-time) and far-field (late-time)
regimes. Each transmitter segment :math:`i` contributes to the field
at the receiver with a weight
.. math::
w_i(t) = \operatorname{erf}\!\left(
\frac{\lambda(t)}{d_i}
\right), \quad
\lambda(t) = \sqrt{\frac{4\,\rho_a(t)\,t}{\mu_0}}
where :math:`d_i` is the perpendicular distance from the receiver to
segment :math:`i`, and :math:`\lambda(t)` is the EM diffusion length.
.. math::
\eta(t) =
\frac{\displaystyle\sum_i w_i^{(r_x)}(t)\,
\Delta H_z^{(i)}(r_x, r_y)}
{\displaystyle\sum_i w_i^{(0)}(t)\,
\Delta H_z^{(i)}(0, 0)}
**Limits**:
* Late time (:math:`\lambda \gg d_i`): :math:`w_i \to 1` for all
segments → :math:`\eta \to \eta_\mathrm{static}` (Biot-Savart limit).
* Early time (:math:`\lambda \ll d_i`): :math:`w_i \propto \lambda/d_i`
→ segments at small :math:`d_i` dominate (near-field limit).
Parameters
----------
rx, ry : float
Receiver offset from loop centre (m).
loop_shape : str
``'square'``, ``'rectangular'``, or ``'circle'``.
loop_dims : tuple
Loop geometry parameters.
t : ndarray (n,)
Gate centre times (s).
rho_a : ndarray (n,)
Current estimate of apparent resistivity (Ω m).
Returns
-------
eta_t : ndarray (n,)
Time-dependent geometry factor.
References
----------
.. [1] Nabighian, M.N. & Macnae, J.C. (1991). SEG, Vol. 2, eq. 4.34.
.. [2] Christiansen, A.V. et al. (2009). *Geophysics*, 74, F35–F46.
"""
from scipy.special import erf
t = np.asarray(t)
n = len(t)
if abs(rx) < 1e-10 and abs(ry) < 1e-10:
return np.ones(n)
# --- get per-segment Hz and distances ---
if loop_shape in ("square", "rectangular"):
if loop_shape == "square":
a = b = float(loop_dims[0]) / 2.0
else:
a = float(loop_dims[0]) / 2.0
b = float(loop_dims[1]) / 2.0
hz_rx, hz_00, dist_rx, dist_00 = _biot_savart_rect_hz_segments(
rx, ry, a, b
)
elif loop_shape == "circle":
r_loop = float(loop_dims[0])
hz_rx, hz_00, dist_rx, dist_00 = _biot_savart_circle_hz_segments(
rx, ry, r_loop
)
else:
a = b = np.sqrt(float(loop_dims[0])) / 2.0
hz_rx, hz_00, dist_rx, dist_00 = _biot_savart_rect_hz_segments(
rx, ry, a, b
)
# --- diffusion length per gate ---
rho_safe = np.maximum(np.asarray(rho_a), 1e-6)
lam = np.sqrt(4.0 * rho_safe * t / MU0) # (n,)
# --- per-segment weights (n_seg, n) ---
d_rx = np.maximum(dist_rx[:, None], 1e-12) # avoid /0
d_00 = np.maximum(dist_00[:, None], 1e-12)
w_rx = erf(lam[None, :] / d_rx) # (n_seg, n)
w_00 = erf(lam[None, :] / d_00)
# --- time-dependent weighted Hz ---
hz_td_rx = (hz_rx[:, None] * w_rx).sum(axis=0) # (n,)
hz_td_00 = (hz_00[:, None] * w_00).sum(axis=0)
with np.errstate(divide="ignore", invalid="ignore"):
eta_t = np.where(
np.abs(hz_td_00) > 1e-30, hz_td_rx / hz_td_00, np.ones(n)
)
return eta_t
def _loop_inner_radius(loop_shape: str, loop_dims: tuple) -> float:
"""
Maximum Tx–Rx offset for the receiver to still lie inside the loop.
* ``'square'`` → half the side length
* ``'rectangular'`` → half the shorter side
* ``'circle'`` → the loop radius
"""
if loop_shape == "circle":
return float(loop_dims[0])
if loop_shape == "square":
return float(loop_dims[0]) / 2.0
# rectangular: min(length, width) / 2
return min(float(loop_dims[0]), float(loop_dims[1])) / 2.0
# ---------------------------------------------------------------------------
# Apparent-resistivity formulas
# ---------------------------------------------------------------------------
def _rho_a_late_time(
dBdt: np.ndarray,
t: np.ndarray,
moment: float,
) -> np.ndarray:
r"""
Late-time apparent resistivity for a coincident / central-loop system.
Derived from the late-time asymptotic response of a magnetic
dipole source over a uniform half-space (Ward & Hohmann 1988,
eq. 4.96; Nabighian & Macnae 1991):
.. math::
\rho_a(t) =
\left(\frac{M\,\mu_0^{5/2}}
{10\,\sqrt{\pi}\,|\partial B_z/\partial t|\,t^{5/2}}
\right)^{2/3}
Parameters
----------
dBdt : ndarray (n,)
:math:`\partial B_z/\partial t` in T s⁻¹.
t : ndarray (n,)
Time gates in seconds.
moment : float
Transmitter magnetic moment :math:`M = I n A` in A m².
Returns
-------
ndarray (n,)
Apparent resistivity in Ω m.
References
----------
.. [1] Ward, S.H. & Hohmann, G.W. (1988). Electromagnetic Theory
for Geophysical Applications. Chapter 4, eq. 4.96.
.. [2] Nabighian, M.N. & Macnae, J.C. (1991). Time domain
electromagnetic prospecting methods. SEG, Vol. 2.
"""
with np.errstate(divide="ignore", invalid="ignore"):
arg = (moment * MU0**2.5) / (
10.0 * np.sqrt(np.pi) * np.abs(dBdt) * t**2.5
)
rho = arg ** (2.0 / 3.0)
return np.where(np.isfinite(rho) & (rho > 0.0), rho, np.nan)
def _rho_a_in_loop(
dBdt: np.ndarray,
t: np.ndarray,
moment: float,
loop_shape: str,
loop_dims: tuple,
rx: float,
ry: float = 0.0,
n_iter: int = 3,
) -> np.ndarray:
r"""
Apparent resistivity for an in-loop receiver using a time-dependent
geometry correction.
Starting from the static Biot-Savart factor
:math:`\eta_0 = H_z^\mathrm{static}(r_x, r_y) / H_z^\mathrm{static}(0,0)`,
the algorithm iteratively refines the per-gate factor
:math:`\eta(t)` via :func:`_in_loop_geometry_factor_td`:
.. math::
\rho_a^{(k+1)}(t) = \left(
\frac{M\,\eta^{(k)}(t)\,\mu_0^{5/2}}
{10\,\sqrt{\pi}\,|\partial B_z/\partial t|\,t^{5/2}}
\right)^{2/3}
At late time :math:`\eta(t) \to \eta_0` (static limit); at early
time the correction accounts for the proximity of the receiver to
individual transmitter segments. Convergence is typically reached
in 2–3 iterations.
Parameters
----------
dBdt : ndarray (n,)
t : ndarray (n,)
moment : float
loop_shape : str
``'square'``, ``'rectangular'``, or ``'circle'``.
loop_dims : tuple
Loop geometry parameters.
rx : float
Receiver x-offset from loop centre (m).
ry : float
Receiver y-offset from loop centre (m). Default 0.
n_iter : int
Number of correction iterations (0 = static only). Default 3.
Returns
-------
ndarray (n,)
Corrected apparent resistivity in Ω m.
References
----------
.. [1] Nabighian & Macnae (1991), Table 4.1.
.. [2] Christiansen, A.V. et al. (2009). *Geophysics*, 74, F35–F46.
"""
# Seed: static correction
eta_0 = _in_loop_geometry_factor(rx, ry, loop_shape, loop_dims)
rho_a = _rho_a_late_time(dBdt, t, moment * eta_0)
if n_iter == 0 or (abs(rx) < 1e-10 and abs(ry) < 1e-10):
return rho_a
for _ in range(n_iter):
# Fill NaN/non-positive ρ_a with median for stable η(t) computation
finite = rho_a[np.isfinite(rho_a) & (rho_a > 0.0)]
fill = float(np.median(finite)) if len(finite) > 0 else 100.0
rho_for_eta = np.where(
np.isfinite(rho_a) & (rho_a > 0.0), rho_a, fill
)
eta_t = _in_loop_geometry_factor_td(
rx, ry, loop_shape, loop_dims, t, rho_for_eta
)
rho_a = _rho_a_late_time(dBdt, t, moment * eta_t)
return rho_a
def _rho_a_offset_loop(
dBdt: np.ndarray,
t: np.ndarray,
moment: float,
offset: float,
) -> np.ndarray:
r"""
Late-time apparent resistivity for a separated-loop system where
the receiver is at horizontal distance ``offset`` from the Tx.
Uses the Ward & Hohmann (1988) late-time formula for a vertical
magnetic dipole measured at horizontal distance :math:`d`:
.. math::
\frac{\partial H_z}{\partial t}
= -\frac{M\,\mu_0^{5/2}\,\sigma^{3/2}}
{20\,\sqrt{\pi}\,d^3\,t^{5/2}}
Solving for apparent resistivity:
.. math::
\rho_a(t) = \left(
\frac{M\,\mu_0^{5/2}}
{20\,\sqrt{\pi}\,|\partial B_z/\partial t|\,d^3\,t^{5/2}}
\right)^{2/3}
Parameters
----------
dBdt : ndarray (n,)
Measured :math:`\partial B_z/\partial t` at the Rx position (T s⁻¹).
t : ndarray (n,)
Time gates in seconds.
moment : float
Transmitter moment :math:`M = I n A` in A m².
offset : float
Horizontal Tx–Rx separation :math:`d` in metres.
Returns
-------
ndarray (n,)
Apparent resistivity in Ω m.
Raises
------
ValueError
If ``offset <= 0``.
Notes
-----
The factor 20 (vs 10 in the central-loop formula) and the
:math:`d^3` denominator arise from the angular geometry of the
off-axis dipole field. The formula is valid when :math:`d` is
large compared to the loop size and large compared to the skin
depth (:math:`d\sqrt{\mu_0 \sigma / (4t)} \ll 1`).
References
----------
.. [1] Ward, S.H. & Hohmann, G.W. (1988). Chapter 4, eq. 4.97.
.. [2] Nabighian, M.N. (1979). Quasi-static transient response of
a conducting half-space. *Geophysics*, 44, 1700–1705.
"""
if offset <= 0.0:
raise ValueError(
f"offset must be > 0 for a separated-loop configuration "
f"(got {offset:.3g} m). Use offset=0 for central-loop."
)
with np.errstate(divide="ignore", invalid="ignore"):
arg = (moment * MU0**2.5) / (
20.0 * np.sqrt(np.pi) * np.abs(dBdt) * (offset**3) * t**2.5
)
rho = arg ** (2.0 / 3.0)
return np.where(np.isfinite(rho) & (rho > 0.0), rho, np.nan)
# ---------------------------------------------------------------------------
# Shared downstream helpers (unchanged from v1)
# ---------------------------------------------------------------------------
def _pseudo_freq(t: np.ndarray, convention: str = "skin_depth") -> np.ndarray:
r"""
Map time gates to equivalent pseudo-frequencies.
``"skin_depth"`` (default) :
:math:`f = 1 / (2\pi t)`
``"diffusion"`` :
:math:`f = 1 / (4\pi t)`
"""
if convention == "skin_depth":
return 1.0 / (2.0 * np.pi * t)
if convention == "diffusion":
return 1.0 / (4.0 * np.pi * t)
raise ValueError(
f"Unknown convention '{convention}'. Use 'skin_depth' or 'diffusion'."
)
def _phase_from_rho(
rho_a: np.ndarray,
freq: np.ndarray,
mode: str = "homogeneous",
) -> np.ndarray:
r"""
Estimate MT phase from the apparent resistivity curve.
``"homogeneous"`` — uniform half-space: :math:`\phi_{xy} = 45°`.
``"weidelt"`` — dispersion relation (Weidelt 1972, Boehl et al.
1983):
.. math::
\phi(\omega) \approx \frac{\pi}{4}
+ \frac{1}{2}\frac{\partial\ln\rho_a}{\partial\ln\omega}
"""
if mode == "homogeneous":
return np.full_like(rho_a, 45.0)
if mode == "weidelt":
log_rho = np.log(rho_a)
log_omega = np.log(2.0 * np.pi * freq)
d_log_rho = np.gradient(log_rho, log_omega)
phase_deg = np.degrees(np.pi / 4.0 + 0.5 * d_log_rho)
phase_deg = np.clip(phase_deg, 0.0, 90.0)
phase_deg[~np.isfinite(rho_a)] = np.nan
return phase_deg
raise ValueError(
f"Unknown phase mode '{mode}'. Use 'homogeneous' or 'weidelt'."
)
def _build_z_array(
rho_a: np.ndarray,
phase_xy_deg: np.ndarray,
freq: np.ndarray,
) -> np.ndarray:
r"""
Construct complex impedance tensor array from apparent resistivity
and phase (1-D earth assumption).
.. math::
|Z_{xy}(\omega)| = \sqrt{\rho_a(\omega)\,\omega\,\mu_0}
"""
omega = 2.0 * np.pi * freq
Z_mag = np.sqrt(rho_a * omega * MU0)
phi_rad = np.radians(phase_xy_deg)
Z_xy = Z_mag * (np.cos(phi_rad) + 1j * np.sin(phi_rad))
n = len(freq)
Z_arr = np.zeros((n, 2, 2), dtype=complex)
Z_arr[:, 0, 1] = Z_xy
Z_arr[:, 1, 0] = -Z_xy
return Z_arr
def _z_error_from_rho_error(
rho_a: np.ndarray,
rho_err: np.ndarray,
freq: np.ndarray,
) -> np.ndarray:
r"""
Propagate apparent-resistivity uncertainty to impedance uncertainty.
.. math::
\delta|Z| = \frac{\delta\rho_a}{2\rho_a} \cdot |Z|
"""
omega = 2.0 * np.pi * freq
Z_mag = np.sqrt(rho_a * omega * MU0)
with np.errstate(divide="ignore", invalid="ignore"):
dZ = np.where(rho_a > 0.0, rho_err / (2.0 * rho_a) * Z_mag, np.nan)
n = len(freq)
Z_err = np.zeros((n, 2, 2))
Z_err[:, 0, 1] = dZ
Z_err[:, 1, 0] = dZ
return Z_err
# Expose build_adjacency from gcn for import convenience; harmless if AI not installed
try:
from ..ai.nets.gcn import build_adjacency
except Exception:
build_adjacency = None # type: ignore[assignment]
# ---------------------------------------------------------------------------
# Waveform deconvolution helpers
# ---------------------------------------------------------------------------
def _waveform_moments(waveform, n_samples: int = 2000) -> tuple[float, float]:
r"""
Compute the effective time shift and gate-rejection window for any
transmitter waveform.
The turn-off impulse response is :math:`w(\tau) = -dI/d\tau \geq 0`.
Two cases are handled:
**Turn-off starts at** :math:`t = 0` (e.g. ``RampWaveform``):
:math:`I(0) > 0`. The normalised first moment is
.. math::
\tau_\mathrm{eff} = \int_0^{\tau_\mathrm{win}} \tau\,w(\tau)\,d\tau
and :math:`\tau_\mathrm{window}` is the 99th-percentile of the
cumulative turn-off. Gates with :math:`t_\mathrm{corr} \leq
\tau_\mathrm{window}` are discarded.
**Turn-off completed before** :math:`t = 0` (e.g. ``HalfSineWaveform``):
:math:`I(0) \approx 0`. The last turn-off in the preceding half-period
is identified and its centroid gives :math:`\tau_\mathrm{eff}` (positive,
measured from centroid to :math:`t = 0`). No gate rejection is applied
(:math:`\tau_\mathrm{window} = 0`).
For a :class:`~pycsamt.tdem.waveform.RampWaveform` the result is
analytic: :math:`\tau_\mathrm{eff} = \tau_r/2`,
:math:`\tau_\mathrm{window} = \tau_r`.
Parameters
----------
waveform :
Any waveform object with a ``current_at(t)`` method and a
``base_frequency`` attribute, or ``None``.
n_samples : int
Number of integration points for numerical waveforms.
Returns
-------
tau_eff : float
Effective time shift (seconds).
tau_window : float
Gate-rejection threshold (seconds).
"""
from .waveform import RampWaveform, SquareWaveform
if waveform is None or isinstance(waveform, SquareWaveform):
return 0.0, 0.0
if isinstance(waveform, RampWaveform):
t_ramp = float(waveform.ramp_off)
return t_ramp / 2.0, t_ramp
hp = 0.5 / float(waveform.base_frequency) # half-period
I0 = float(waveform.current_at(np.array([0.0]))[0])
if I0 > 0.01:
# Turn-off begins at (or just after) t=0 — sample [0, 2×hp]
tau_grid = np.linspace(0.0, 2.0 * hp, n_samples + 1)
I = waveform.current_at(tau_grid)
dI = np.diff(I)
tau_mid = 0.5 * (tau_grid[:-1] + tau_grid[1:])
w = np.where(dI < 0.0, -dI, 0.0)
total = float(w.sum())
if total <= 0.0:
return 0.0, 0.0
w_norm = w / total
tau_eff = float(np.dot(tau_mid, w_norm))
cumulative = np.cumsum(w_norm)
idx = int(np.searchsorted(cumulative, 0.99))
tau_window = float(tau_mid[min(idx, len(tau_mid) - 1)])
return tau_eff, tau_window
else:
# Turn-off completed before t=0 — find the last turn-off in [-2hp, 0]
tau_grid = np.linspace(-2.0 * hp, 0.0, n_samples + 1)
I = waveform.current_at(tau_grid)
dI = np.diff(I)
tau_mid = 0.5 * (tau_grid[:-1] + tau_grid[1:])
w = np.where(dI < 0.0, -dI, 0.0)
# Isolate the last contiguous turn-off (nearest to t=0)
inc_idx = np.where(dI > 0.0)[0]
start = int(inc_idx[-1]) + 1 if len(inc_idx) > 0 else 0
w_last = np.zeros_like(w)
w_last[start:] = w[start:]
total = float(w_last.sum())
if total <= 0.0:
return 0.0, 0.0
w_norm = w_last / total
# centroid is at some tau < 0; tau_eff is the positive lag to t=0
tau_eff = -float(np.dot(tau_mid, w_norm))
return tau_eff, 0.0 # turn-off already complete: no gate rejection
def _apply_waveform_correction(
dBdt: np.ndarray,
t: np.ndarray,
waveform,
error: np.ndarray | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]:
r"""
First-order waveform deconvolution for any transmitter waveform.
For a transmitter with a finite turn-off, the measured :math:`\partial
B_z/\partial t` is the convolution of the ideal step-off response
:math:`b(t)` with the normalised turn-off impulse
:math:`w(\tau)=-dI/d\tau`:
.. math::
V_\mathrm{meas}(t) = \int_0^{\tau_\mathrm{max}} w(\tau)\,b(t-\tau)\,d\tau
A first-order Taylor expansion (Fitterman & Stewart 1986) yields
.. math::
b(t) \approx V_\mathrm{meas}(t + \tau_\mathrm{eff})
\left(1 + \frac{\tau_\mathrm{eff}}{2\,t}\right)
where :math:`\tau_\mathrm{eff}` is the first moment of :math:`w`.
Gates measured within the turn-off window (:math:`t_\mathrm{corr}
\leq \tau_\mathrm{window}`) are discarded.
Waveform-type dispatch:
* :class:`~pycsamt.tdem.waveform.SquareWaveform` — no correction (ideal).
* :class:`~pycsamt.tdem.waveform.RampWaveform` — analytic
:math:`\tau_\mathrm{eff}=\tau_r/2`, recovering the classic
Fitterman & Stewart formula exactly.
* :class:`~pycsamt.tdem.waveform.HalfSineWaveform` and
:class:`~pycsamt.tdem.waveform.CustomWaveform` — :math:`\tau_\mathrm{eff}`
computed numerically via :func:`_waveform_moments`.
Parameters
----------
dBdt : np.ndarray
Measured :math:`\partial B_z/\partial t` gate values.
t : np.ndarray
Gate centre times in seconds.
waveform :
Waveform descriptor or ``None``.
error : np.ndarray or None
Per-gate measurement uncertainty; corrected identically to ``dBdt``.
Returns
-------
dBdt_corr : np.ndarray
t_corr : np.ndarray
error_corr : np.ndarray or None
"""
tau_eff, tau_window = _waveform_moments(waveform)
if tau_eff <= 0.0:
return dBdt, t, error
t_corr = t + tau_eff
amp_corr = 1.0 + tau_eff / (2.0 * t_corr)
valid = t_corr > tau_window
err_corr = error[valid] if error is not None else None
return dBdt[valid] * amp_corr[valid], t_corr[valid], err_corr
# ---------------------------------------------------------------------------
# LateTimeTransform
# ---------------------------------------------------------------------------
_CONFIG_CENTRAL = "central"
_CONFIG_IN_LOOP = "in_loop"
_CONFIG_OFFSET = "offset"
# ---------------------------------------------------------------------------
# Fourier / Kramers-Kronig numerical helpers
# ---------------------------------------------------------------------------
def _cosine_transform_1d(
g: np.ndarray,
t: np.ndarray,
omega_arr: np.ndarray,
n_interp: int = 0,
) -> np.ndarray:
r"""
Half-range Fourier cosine transform via trapezoidal quadrature in
log-time space.
.. math::
F_c[g](\omega) = \int_0^\infty g(t)\,\cos(\omega t)\,\mathrm{d}t
\approx \int_{\ln t_0}^{\ln t_N}
g(t)\,\cos(\omega t)\,t\,\mathrm{d}(\ln t)
Parameters
----------
g : ndarray (n,)
Integrand at discrete time gates.
t : ndarray (n,)
Strictly positive, monotone increasing time gates in seconds.
omega_arr : ndarray (m,)
Angular frequencies [rad/s] at which to evaluate.
n_interp : int
If > ``len(t)``, interpolate *g* onto a denser log-spaced grid
before integration. Default 0 (no interpolation).
Returns
-------
ndarray (m,)
"""
if n_interp > len(t):
t_fine = np.exp(
np.linspace(np.log(t[0]), np.log(t[-1]), int(n_interp))
)
log_abs_g = np.log(np.maximum(np.abs(g), 1e-300))
g_fine = np.sign(np.interp(np.log(t_fine), np.log(t), g)) * np.exp(
np.interp(np.log(t_fine), np.log(t), log_abs_g)
)
else:
t_fine = t
g_fine = g
log_t = np.log(t_fine)
result = np.empty(len(omega_arr))
for i, w in enumerate(omega_arr):
integrand = g_fine * np.cos(w * t_fine) * t_fine # ×t : d(log t)→dt
result[i] = np.trapezoid(integrand, log_t)
return result
def _kramers_kronig_re(
omega: np.ndarray,
im_k: np.ndarray,
) -> np.ndarray:
r"""
Compute :math:`\mathrm{Re}[K(\omega)]` from
:math:`\mathrm{Im}[K(\omega)]` via the one-sided K-K relation:
.. math::
\mathrm{Re}[K(\omega)] =
\frac{2}{\pi}\,\mathrm{P.V.}
\int_0^\infty
\frac{\omega'\,\mathrm{Im}[K(\omega')]}{\omega'^2 - \omega^2}
\,\mathrm{d}\omega'
The Cauchy principal value is handled by interpolating across the
singular point :math:`\omega' = \omega` in log-space.
Parameters
----------
omega : ndarray (n,) angular frequencies, log-spaced recommended
im_k : ndarray (n,) imaginary part of the kernel
Returns
-------
ndarray (n,) real part of the kernel
"""
n = len(omega)
log_omega = np.log(omega)
num_base = omega**2 * im_k # ω'^2 Im[K] — numerator in log-ω integrand
re_k = np.zeros(n)
for i in range(n):
w2 = omega[i] ** 2
den = omega**2 - w2 # zero only at j == i
tol = 1e-10 * max(w2, 1e-100)
mask_ok = np.abs(den) > tol
integrand = np.where(
mask_ok, num_base / np.where(mask_ok, den, 1.0), 0.0
)
bad = np.where(~mask_ok)[0]
for j in bad:
left = j - 1 if j > 0 else j + 1
right = j + 1 if j < n - 1 else j - 1
integrand[j] = 0.5 * (integrand[left] + integrand[right])
re_k[i] = (2.0 / np.pi) * np.trapezoid(integrand, log_omega)
return re_k
# ---------------------------------------------------------------------------
# FourierTransform
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# TEMtoEDI — high-level dispatcher
# ---------------------------------------------------------------------------
[docs]
class TEMtoEDI(PyCSAMTObject):
r"""
Convert one or more TEM soundings to a
:class:`~pycsamt.seg.collection.EDICollection`.
Wraps either :class:`LateTimeTransform` (default) or
:class:`FourierTransform` and writes one synthetic EDI file per
sounding into ``out_dir``.
All three loop configurations (central, in-loop, offset) are handled
automatically via :class:`LateTimeTransform`.
Parameters
----------
method : str
``"late_time"`` (default) or ``"fourier"``.
freq_convention : str
``"skin_depth"`` (default) or ``"diffusion"``.
phase_mode : str
``"homogeneous"`` (default) or ``"weidelt"``.
loop_geometry_correction : bool
Forward to :class:`LateTimeTransform`. Default ``True``.
out_dir : str or Path
EDI output directory. Default ``"edi_out/tem"``.
verbose : int
Verbosity. Default 0.
Examples
--------
Central-loop:
>>> from pycsamt.tdem import TEMSounding, TEMtoEDI
>>> import numpy as np
>>> t = np.logspace(-5, -2, 25)
>>> dBdt = 5e-5 * t ** (-5.0 / 2.0)
>>> snd = TEMSounding(t, dBdt, current=8.0, tx_area=1e4,
... station_name="S01")
>>> conv = TEMtoEDI(method="late_time")
>>> coll = conv.transform(snd)
Offset (separated) loop — just set ``offset``:
>>> snd_off = TEMSounding(t, dBdt, current=8.0, tx_area=1e4,
... offset=500.0, loop_dims=(100.0,),
... station_name="S01")
>>> coll_off = conv.transform(snd_off)
"""
def __init__(
self,
method: str = "late_time",
freq_convention: str = "skin_depth",
phase_mode: str = "homogeneous",
loop_geometry_correction: bool = True,
out_dir: str | Path = "edi_out/tem",
verbose: int = 0,
logger: object | None = None,
) -> None:
self.method = method.lower()
self.out_dir = out_dir
self.verbose = int(verbose)
self.logger = logger
if self.method == "late_time":
self._transformer = LateTimeTransform(
freq_convention=freq_convention,
phase_mode=phase_mode,
loop_geometry_correction=loop_geometry_correction,
verbose=verbose,
logger=logger,
)
elif self.method == "fourier":
self._transformer = FourierTransform(
verbose=verbose,
logger=logger,
)
else:
raise ValueError(
f"Unknown method '{method}'. Use 'late_time' or 'fourier'."
)
[docs]
def save(self, sounding_or_soundings, path=None):
"""Transform and write EDI files. Returns list of written paths."""
from pathlib import Path as _Path
out = _Path(self.out_dir if path is None else path)
out.mkdir(parents=True, exist_ok=True)
soundings = (
[sounding_or_soundings]
if isinstance(sounding_or_soundings, TEMSounding)
else list(sounding_or_soundings)
)
written = []
for snd in soundings:
result = self._transformer.transform(snd)
edi = self._result_to_edifile(result)
name = (result["station_name"] or "tem_site") + ".edi"
fp = str(out / name)
edi.write(new_edifn=name, savepath=str(out))
written.append(fp)
if self.verbose:
cfg = result.get("loop_config", "?")
print(
f" Wrote {fp} ({result['freq'].size} freq, cfg={cfg})"
)
return written
# ------------------------------------------------------------------
# Internal helpers
# ------------------------------------------------------------------
def _result_to_edifile(self, result: dict):
from ..seg.edi import EDIFile
edi = EDIFile()
edi.Z.freq = result["freq"]
edi.Z.z = result["Z"]
z_err = self._finite_z_error(result.get("Z_err"))
if z_err is not None:
edi.Z.z_err = z_err
self._set_head(edi, result)
return edi
@staticmethod
def _finite_z_error(z_err):
"""Return finite non-negative impedance errors when present."""
if z_err is None:
return None
arr = np.asarray(z_err, dtype=float)
if not np.isfinite(arr).any():
return None
arr = np.where(np.isfinite(arr), arr, 0.0)
arr = np.maximum(arr, 0.0)
return arr
@staticmethod
def _set_head(edi, result: dict) -> None:
cfg = result.get("loop_config", "central")
try:
from ..seg.heads import Head
head = Head()
head.dataid = result["station_name"] or "TEM_SITE"
head.lon = float(result.get("x", 0.0))
head.lat = float(result.get("y", 0.0))
head.elev = float(result.get("elevation", 0.0))
head.acqby = "pyCSAMT.tdem"
head.fileby = "pyCSAMT.tdem"
edi.add_section("head", head)
except Exception:
pass
try:
from ..seg.heads import Info
info = Info()
info.info_list = [
"Data origin: TEM (time-domain EM) transformed to MT equivalent",
"Transform: pycsamt.tdem.transform.TEMtoEDI",
f"Loop configuration: {cfg}",
"Z_xx = Z_yy = 0 (1-D assumption); Z_yx = -Z_xy",
]
edi.add_section("info", info)
except Exception:
pass
@staticmethod
def _make_collection(edis: list):
from ..seg.collection import EDICollection
return EDICollection(items=edis)
