From 8ce1c87954cd0ea7cfaaa8b7705fe9aa675640d1 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Fri, 23 Jan 2026 19:19:24 +0100 Subject: [PATCH 01/14] - Implement waveform through signal - Implement bandpass through user func - Update examples: - IP and VRM - TEM-FAST - WalkTEMT --- CHANGELOG.rst | 3 + empymod/model.py | 33 ++- empymod/utils.py | 87 +++++- examples/time_domain/ip_vrm.py | 428 ++++++++-------------------- examples/time_domain/tem_temfast.py | 378 ++++++------------------ examples/time_domain/tem_walktem.py | 391 ++++++++----------------- 6 files changed, 439 insertions(+), 881 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 4c6123a..f239feb 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -15,6 +15,9 @@ v2.5.x - Modelling routines: + - Arbitrary waveforms: Time-domain modelling can new be done not only for + impulse (``signal=0``), step-on (``signal=1``), and step-off + (``signal=-1``), but for any arbitrary piecewise linear waveform. - Merged ``loop`` into ``bipole``; there is no need for a special routine. Simply use ``bipole`` with ``msrc='b'`` and ``mrec=True`` for the same effect. As a result, ``empymod.utils.check_bipole`` changed its signature. diff --git a/empymod/model.py b/empymod/model.py index 7db74f7..14c20d5 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -325,6 +325,9 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, If True, the output is squeezed. If False, the output will always be of ``ndim=3``, (nfreqtime, nrec, nsrc). + bandbass : dict or None, default: None + TODO + Returns ------- @@ -404,12 +407,16 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, Out[2]: (1.6880934577857306e-10-3.083031298956568e-10j) """ + # Get kwargs with defaults. out = get_kwargs( - ['verb', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'squeeze'], - [2, 'dlf', {}, 'dlf', {}, False, None, True], kwargs, + [ + 'verb', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'squeeze', + 'bandpass', + ], + [2, 'dlf', {}, 'dlf', {}, False, None, True, None], kwargs, ) - verb, ht, htarg, ft, ftarg, xdirect, loop, squeeze = out + verb, ht, htarg, ft, ftarg, xdirect, loop, squeeze, bandpass = out # === 1. LET'S START ============ t0 = printstartfinish(verb) @@ -420,7 +427,8 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, if signal is None: freq = freqtime else: - time, freq, ft, ftarg = check_time(freqtime, signal, ft, ftarg, verb) + outtime = check_time(freqtime, signal, ft, ftarg, verb, True) + time, freq, ft, ftarg, signal, waveform = outtime # Check layer parameters model = check_model(depth, res, aniso, epermH, epermV, mpermH, mpermV, @@ -598,6 +606,10 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, # In case of QWE/QUAD, print Warning if not converged conv_warning(conv, htarg, 'Hankel', verb) + # Apply user-provided bandpass filter + if isinstance(bandpass, dict) and 'func' in bandpass: + bandpass['func'](bandpass, locals()) + # Do f->t transform if required if signal is not None: EM, conv = tem(EM, EM[0, :], freq, time, signal, ft, ftarg) @@ -605,8 +617,21 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, # In case of QWE/QUAD, print Warning if not converged conv_warning(conv, ftarg, 'Fourier', verb) + # Apply waveform + if waveform: + map_time = waveform['map_time'][:, :, :] + wave_weight = waveform['wave_weight'][:, :, None] + wave_index = waveform['wave_index'][:, :, None] + gauss_weight = waveform['gauss_weight'][None, :, None, None] + + # Convert to time domain and apply waveform + EM = np.sum(np.sum(EM[map_time] * gauss_weight, axis=1) * + wave_weight * wave_index, axis=1) + # Reshape for number of sources EM = EM.reshape((-1, nrec, nsrc), order='F') + + # Squeeze empty dimensions if squeeze: EM = np.squeeze(EM) diff --git a/empymod/utils.py b/empymod/utils.py index b700888..cfcd6e1 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -56,7 +56,7 @@ 'check_bipole', 'check_ab', 'check_solution', 'get_abs', 'get_geo_fact', 'get_azm_dip', 'get_off_ang', 'get_layer_nr', 'printstartfinish', 'conv_warning', 'set_minimum', 'get_minimum', - 'Report'] + 'Report', 'check_waveform'] # 0. General settings @@ -961,7 +961,7 @@ def check_loop(loop, ht, htarg, verb): return loop_freq, loop_off -def check_time(time, signal, ft, ftarg, verb): +def check_time(time, signal, ft, ftarg, verb, new=False): r"""Check time domain specific input parameters. This check-function is called from one of the modelling routines in @@ -1005,6 +1005,11 @@ def check_time(time, signal, ft, ftarg, verb): """ + if new and isinstance(signal, dict): + time, signal, waveform = check_waveform(time, **signal) + else: + waveform = None + # Check time and input signal time = check_time_only(time, signal, verb) @@ -1249,7 +1254,11 @@ def check_time(time, signal, ft, ftarg, verb): if args and verb > 0: print(f"* WARNING :: Unknown ftarg {args} for method '{ft}'") - return time, freq, ft, targ + # Make backwards compatible with old signature possibility? + if new: + return time, freq, ft, targ, signal, waveform + else: + return time, freq, ft, targ def check_time_only(time, signal, verb): @@ -1301,6 +1310,70 @@ def check_time_only(time, signal, verb): return time +def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): + + # The waveform function (here and in model.bipole) is based on work + # from Kerry Key. + + # Waveform segments. + dt = np.diff(nodes) + dI = np.diff(amplitudes) + dIdt = dI/dt + + # Gauss-Legendre Quadrature; 3 is generally good enough. + # (Roots/weights could be cached.) + g_x, gauss_weight = sp.special.roots_legendre(nquad) + + # Pre-allocate + comp_time = np.zeros((time.size, nquad, dIdt.size)) + wave_weight = np.zeros((time.size, dIdt.size)) + wave_index = np.zeros((time.size, dIdt.size), dtype=bool) + + # Loop over wave segments. + wi = 0 + for i, cdIdt in enumerate(dIdt): + + # We only have to consider segments with a change of current. + if cdIdt == 0.0: + continue + + # If wanted time is before a wave element, ignore it. + wave_ind = nodes[i] < time + if wave_ind.sum() == 0: + continue + + # Start and end for this wave-segment for all times. + ta = time[wave_ind] - nodes[i] + tb = time[wave_ind] - nodes[i+1] + + # If wanted time is within a wave element, we cut the element. + tb[nodes[i+1] > time[wave_ind]] = 0.0 # Cut elements + + # Gauss-Legendre for this wave segment. See + # https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval + # for the change of interval, which makes this a bit more complex. + comp_time[:, :, wi] = np.outer((tb-ta)/2, g_x)+(ta+tb)[:, None]/2 + + wave_weight[:, wi] = (tb-ta)/2*cdIdt + wave_index[:, wi] = wave_ind + + wi += 1 + + comp_time = comp_time[:, :, :wi] + wave_weight = wave_weight[:, :wi] + wave_index = wave_index[:, :wi] + + comp_time_flat, map_time = np.unique(comp_time, return_inverse=True) + waveform = { + 'map_time': map_time, + 'wave_weight': wave_weight, + 'wave_index': wave_index, + 'gauss_weight': gauss_weight, + } + + return comp_time_flat, signal, waveform + + def check_solution(solution, signal, ab, msrc, mrec): r"""Check required solution with parameters. @@ -1820,7 +1893,7 @@ def get_kwargs(names, defaults, kwargs): - ALL functions: src, rec, res, aniso, epermH, epermV, mpermH, mpermV, verb - ONLY gpr: cf, gain - - ONLY bipole: msrc, srcpts + - ONLY bipole: msrc, srcpts, bandpass - ONLY dipole_k: freq, wavenumber - ONLY analytical: solution - ONLY bipole, loop: mrec, recpts, strength @@ -1849,9 +1922,9 @@ def get_kwargs(names, defaults, kwargs): """ # Known keys (excludes keys present in ALL routines). known_keys = set([ - 'depth', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'signal', - 'ab', 'freqtime', 'freq', 'wavenumber', 'solution', 'cf', 'gain', - 'msrc', 'srcpts', 'mrec', 'recpts', 'strength', 'squeeze' + 'depth', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'signal', + 'ab', 'freqtime', 'freq', 'wavenumber', 'solution', 'cf', 'gain', + 'msrc', 'srcpts', 'mrec', 'recpts', 'strength', 'squeeze', 'bandpass', ]) # Loop over wanted parameters. diff --git a/examples/time_domain/ip_vrm.py b/examples/time_domain/ip_vrm.py index f6b0ff8..85c9c47 100644 --- a/examples/time_domain/ip_vrm.py +++ b/examples/time_domain/ip_vrm.py @@ -10,11 +10,10 @@ """ import empymod import numpy as np -import scipy as sp import matplotlib.pyplot as plt from scipy.constants import mu_0 -plt.style.use('ggplot') +plt.style.use("ggplot") ############################################################################### @@ -24,260 +23,101 @@ # Loops # ''''' # -# Create a square loop source of 400x400 m, and two Z-component receivers, one -# outside and one inside the loop; all at the surface (z=0). - -# Create dipoles: [x0, x1, y0, y1, z0, z1] -src_x = np.r_[ - np.zeros(10), np.arange(10), np.ones(10)*10, np.arange(10, -1, -1) -]*40 - 200 -src_y = np.r_[ - np.arange(10), np.ones(10)*10, np.arange(10, -1, -1), np.zeros(10) -]*40 - 200 -src_dipole = [src_x[:-1], src_x[1:], src_y[:-1], src_y[1:], 0, 0] - -# Receiver locations: One outside, one inside; vertical +# Create a square loop source of 400 x 400 m, and two Z-component receivers, +# one outside and one inside the loop; all at the surface (z=0). +# +# **Note**: Take care of the direction of the loop; defining it +# counterclockwise will yield responses that are opposite (factor -1) to a loop +# defined clockwise (following Farady's law of induction). + +# Src: x0, x1, y0, y1, z0, z1 +src_x = np.array([-200, -200, 200, 200, -200]) +src_y = np.array([-200, 200, 200, -200, -200]) +src_dipole = [src_x[:-1], src_x[1:], src_y[:-1], src_y[1:], 0, 0] +# Rec: x, y, z, azm, dip rec = [[-400., 0], [0, 0], [0, 0], 0, 90] # Plot the loop fig, ax = plt.subplots(constrained_layout=True) # Source loop -ax.plot(src_x[::10], src_y[::10], 'ko', ms=10, label='Loop corners') -ax.plot(src_x, src_y, 'C2.-', lw=2, label='LineCurrent path') +ax.plot(src_x[:-1], src_y[:-1], "ko", ms=10, label="Loop corners") +ax.plot(src_x, src_y, "C2.-", lw=2, label="Loop path") # Receiver locations -ax.plot(rec[0][0], rec[1][0], 's', label='Outside Rx') -ax.plot(rec[0][1], rec[1][1], 's', label='Inside Rx') +ax.plot(rec[0][0], rec[1][0], "s", label="Outside Rx") +ax.plot(rec[0][1], rec[1][1], "s", label="Inside Rx") -ax.set_xlabel('X (m)') -ax.set_ylabel('Y (m)') -ax.set_title('Survey layout') +ax.set_xlabel("X (m)") +ax.set_ylabel("Y (m)") +ax.set_title("Survey layout") ax.legend() -ax.set_aspect('equal') +ax.set_aspect("equal") ############################################################################### # Trapezoid Waveform # '''''''''''''''''' -def current(times, nodes): - """Small helper routine to get the waveform current for the given times.""" - # Off by default - out = np.zeros(times.size) - - # Ramp on - i = (times >= nodes[0]) * (times <= nodes[1]) - out[i] = (1.0 / (nodes[1] - nodes[0])) * (times[i] - nodes[0]) - - # On - i = (times > nodes[1]) * (times < nodes[2]) - out[i] = 1 - - # Ramp off - i = (times >= nodes[2]) * (times <= nodes[3]) - out[i] = 1 - (1.0 / (nodes[3] - nodes[2])) * (times[i] - nodes[2]) - - return out - - # On-time negative, t=0 at end of ramp-off # Quarter period for 50 % duty cycle source_frequency_hz = 0.25 on_time_s = 1 / (source_frequency_hz * 4) # Time channels: off-time only, positive times -# 25 channels from 0.1 ms to 1000 ms for 0.25 Hz 50% duty cycle +# 25 channels from 0.1 ms to 1000 ms for 0.25 Hz 50 % duty cycle times = np.logspace(-4, np.log10(on_time_s), 25) # Waveform -nodes_times = np.array([-1, -0.999, -0.001, 0]) -nodes_current = np.array([0., 1, 1, 0]) -waveform_times = np.linspace(nodes_times[0] - 0.1, times[-1] + 0.1, 100000) -waveform_current = current(waveform_times, nodes_times) - +nodes = np.array([-1, -0.999, -0.001, 0]) +amplitudes = np.array([0., 1, 1, 0]) -############################################################################### -print("Waveform details:") -print( - f" Time channels: {len(times)} channels from {times[0]*1000:.2f} ms to" - f"{times[-1]*1000:.2f} ms (all off-time)" -) print( - f" Waveform: on-time from {nodes_times[0]:.3f}s to" - f"{nodes_times[3]:.3f}s" + "Waveform details:\n" + f" Time channels: {len(times)} channels " + f"from {times[0]*1000:.2f} ms to {times[-1]:.2f} s (all off-time)\n" + " Waveform:\n" + f" On-time from {nodes[0]:.1f} s to {nodes[3]:.1f} s\n" + f" Ramp on: {nodes[0]*1e3:.1f} to {nodes[1]*1e3:.1f} ms\n" + f" Ramp off: {nodes[2]*1e3:.1f} to {nodes[3]*1e3:.1f} ms\n" ) -print(f" Ramp on: {nodes_times[0]*1e3:.3f} to {nodes_times[1]*1e3:.3f} ms") -print(f" Ramp off: {nodes_times[2]*1e3:.3f} to {nodes_times[3]*1e3:.3f} ms") -print(" Off time: 0.000 ms") - -# Plot waveform and time channels -fig, axs = plt.subplots( - 1, 2, figsize=(14, 4), sharey=True, constrained_layout=True) - -for ax in axs: - # Plot waveform - ax.plot(waveform_times * 1e3, waveform_current, 'b-', linewidth=2, - label='Waveform') - ax.axhline(0, color='k', linestyle='--', linewidth=0.5) - ax.axvline(0, color='k', linestyle='--', linewidth=0.5, - label='t=0 (end of ramp-off)') - - # Mark waveform time nodes - ax.plot(nodes_times * 1e3, nodes_current, 'ro', markersize=8, - label='Waveform nodes', zorder=5) - # Mark time channels - ax.plot(times * 1000, np.zeros_like(times), 'g|', markersize=10, - label='Time channels', zorder=10) - # Formatting - ax.set_xlabel('Time (ms)') +############################################################################### +# Plot waveform and time channels +fig, ax = plt.subplots(1, 1, figsize=(8, 5), constrained_layout=True) -fig.suptitle('Waveform and Time Channels') -axs[0].set_ylabel('Normalized Current') -axs[0].legend() -axs[1].set_xscale('symlog', linthresh=0.4, linscale=0.5) +# Waveform +ax.plot(np.r_[-1.5, nodes, 1.5], np.r_[0, amplitudes, 0], + "C0-", lw=2, label="Waveform") +# Nodes +ax.plot(nodes, amplitudes, "C1o", markersize=8, label="Waveform nodes") + +# Mark time channels +ax.plot(times, np.zeros(times.size), "k|", ms=10, label="Time channels") + +# Formatting +ax.set_xlabel("Time (s)") +ax.set_title("Waveform and Time Channels") +ax.set_ylabel("Normalized Current") +ax.set_xlim([-1.3, 1.3]) +ax.legend() ############################################################################### -# Waveform Functions -# ------------------ +# Viscous Remanent Magnetization (VRM) +# ------------------------------------ # -# These functions handle the trapezoid waveform convolution for the -# simulations. They are adapted from the `WalkTEM example -# `_. +# This function implements the following VRM model # -# Key differences between B field and dB/dt: +# .. math:: # -# - For B field: dipole source strength is `mu_0 * loop_current` -# - For dB/dt: dipole source strength is `loop_current`, and we multiply by -# `i*omega*mu_0` before frequency-to-time conversion - -def get_time(time_channels, nodes_times): - """ - Compute required times for waveform convolution. - - Because of the arbitrary waveform, we need to compute some times before and - after the actually wanted times for interpolation of the waveform. - - time_req : ndarray - Required times for computation (incl. extra points for interpolation) - """ - t_log = np.log10(time_channels) - # Add a point at the minimum time channel minus the time step, but don't go - # lower than t=0 (end of ramp) - tmin = np.max([t_log[0] - (t_log[1] - t_log[0]), -10]) - # Add a point at the maximum time channel plus the time step - tmax = t_log[-1] + (t_log[-1] - t_log[-2]) - return np.logspace(tmin, tmax, time_channels.size + 2) - - -############################################################################### - -def apply_waveform_to_signal( - times, resp, time_channels, wave_times, wave_amp, nquad=3): - """ - Apply a source waveform to the signal. - - Modified from empymod WalkTEM example. - """ - # Interpolate on log. - PP = sp.interpolate.InterpolatedUnivariateSpline(np.log10(times), resp) - - # Wave time steps. - dt = np.diff(wave_times) - dI = np.diff(wave_amp) - dIdt = dI / dt - - # Gauss-Legendre Quadrature; 3 is generally good enough. - g_x, g_w = sp.special.roots_legendre(nquad) - - # Pre-allocate output. - resp_wanted = np.zeros_like(time_channels) - - # Loop over wave segments. - for i, cdIdt in enumerate(dIdt): - # We only have to consider segments with a change of current. - if cdIdt == 0.0: - continue - - # If wanted time is before a wave element, ignore it. - ind_a = wave_times[i] < time_channels - if ind_a.sum() == 0: - continue - - # If wanted time is within a wave element, we cut the element. - ind_b = wave_times[i + 1] > time_channels[ind_a] - - # Start and end for this wave-segment for all times. - ta = time_channels[ind_a] - wave_times[i] - tb = time_channels[ind_a] - wave_times[i + 1] - tb[ind_b] = 0.0 # Cut elements - - # Gauss-Legendre for this wave segment. - logt = np.log10(np.outer((tb - ta) / 2, g_x) + (ta + tb)[:, None] / 2) - fact = (tb - ta) / 2 * cdIdt - resp_wanted[ind_a] += fact * np.sum(np.array(PP(logt) * g_w), axis=1) - - return resp_wanted - - -############################################################################### - -def convert_freq_to_time(EM, freq, time, ft, ftarg, time_channels, nodes_times, - waveform_current, compute_B_field=True): - """ - Convert frequency-domain response to time domain and apply waveform. - - Parameters - ---------- - EM : ndarray - Frequency-domain EM response - freq : ndarray - Frequencies - time : ndarray - Time array for conversion - ft : str - Transform type - ftarg : dict - Transform arguments - time_channels : ndarray - Desired output times - nodes_times : ndarray - Waveform time nodes - waveform_current : ndarray - Waveform current values at time nodes - compute_B_field : bool - If True, compute B field (T). If False, compute dB/dt (T/s). - - Returns - ------- - resp_wanted : ndarray - Time-domain response at time_channels - """ - if not compute_B_field: - # For dB/dt: multiply by i*omega*mu_0 to convert H to dB/dt in - # frequency domain - EM *= 2j * np.pi * freq * mu_0 - - # Convert to time domain - delay_rst = 0 - EM, _ = empymod.model.tem(EM[:, None], np.array([1]), freq, time + - delay_rst, 1, ft, ftarg) - EM = np.squeeze(EM) - - # Apply waveform - return apply_waveform_to_signal(time, EM, time_channels, nodes_times, - waveform_current) - - -############################################################################### -# Viscous Remanent Magnetization (VRM) Function -# --------------------------------------------- +# \chi (\omega ) = \chi + \Delta \chi \left[ +# 1 - \frac{1}{\log (\tau_2 / \tau_1 )} +# \log \left( \frac{1 + i\omega \tau_2}{1 + i\omega \tau_1} \right) +# \right] \ , \qquad\qquad\qquad (1) # -# This function implements VRM modeling, which is supported but not computed -# within ``empymod``. +# where :math:`\chi` is the viscous susceptibility. def vrm_from_mu(inp, p_dict): """ @@ -312,13 +152,21 @@ def vrm_from_mu(inp, p_dict): ############################################################################### -# Cole-Cole Function -# ------------------ +# Pelton Cole-Cole +# ---------------- +# +# This function implements the following Pelton Cole-Cole IP model for the +# resistivities +# +# .. math:: +# +# \rho(\omega) = \rho_\infty \left[1 + \frac{m}{(1 - m)(1 + +# (i\omega\tau)^C)} \right]\ , \qquad\qquad\qquad (2) # -# This function implements Cole-Cole IP modeling, which is supported but not -# computed within ``empymod``. For more info on the Pelton model refer to the -# `IP example -# `_. +# where :math:`m` is the intrinsic chargeablitiy. For more information on the +# Cole-Cole model refer to the IP example +# :ref:`sphx_glr_gallery_tdomain_cole_cole_ip.py`; note that this model is +# slightly different from the one presented there. def pelton_cole_cole_model(inp, p_dict): """ @@ -334,8 +182,9 @@ def pelton_cole_cole_model(inp, p_dict): """ # Compute complex resistivity from Pelton et al. iotc = np.outer(2j * np.pi * p_dict["freq"], inp["tau"]) ** inp["c"] + # Version using equation 16 of Tarasov & Titov (2013) - rhoH = inp["res"] * (1 + (inp["m"] / (1 - inp["m"])) * (1 / (1 + iotc))) + rhoH = inp["res"] * (1 + inp["m"] / (1 - inp["m"]) / (1 + iotc)) rhoV = rhoH * p_dict["aniso"] ** 2 # Add electric permittivity contribution @@ -345,96 +194,41 @@ def pelton_cole_cole_model(inp, p_dict): return etaH, etaV -############################################################################### -# Custom empymod routine for IP and VRM -# ------------------------------------- - -def simulate_empymod(inp, res, times, compute_B_field=True, loop_current=1.0, - apply_vrm=False, apply_cole_cole=False): - """ - Simulate TDEM response using empymod with trapezoid waveform. Supports VRM - and Cole-Cole IP modeling. - """ - # Get required times for computation - time = get_time(times, nodes_times) - - # Get required frequencies - time, freq, ft, ftarg = empymod.utils.check_time( - time=time, - signal=-1, # Switch-off response - ft="dlf", - ftarg={"dlf": "key_81_2009"}, - verb=inp.get('verb', 0), - ) - - # Source strength depends on output field type - # For B field: mu_0 * current returns B in frequency domain - # For dB/dt: just current, then multiply by i*omega*mu_0 later - strength = mu_0 * loop_current if compute_B_field else loop_current - - # Build empymod model dict with VRM and/or Cole-Cole hooks as needed - res_dict = {"res": res} - - if apply_vrm: - res_dict = {**res_dict, "func_zeta": vrm_from_mu, **apply_vrm} - - if apply_cole_cole: - res_dict = {**res_dict, "func_eta": pelton_cole_cole_model, - **apply_cole_cole} - - # Compute frequency-domain response (summed over source elements) - EM_loop = empymod.model.bipole( - **inp, - res=res_dict, - freqtime=freq, - signal=None, - mrec=True, - strength=strength, - epermH=np.r_[0.0, np.ones(len(res)-1)], - msrc=False, - srcpts=3, - htarg={"dlf": "key_101_2009", "pts_per_dec": -1}, - ).sum(axis=-1) - - # Convert to time domain and apply waveform - nrec = EM_loop.shape[1] - responses = np.zeros((nrec, times.size)) - - for i_rx in range(nrec): - responses[i_rx, :] = convert_freq_to_time( - EM_loop[:, i_rx], freq, time, ft, ftarg, times, nodes_times, - nodes_current, compute_B_field=compute_B_field - ) - - return responses - - ############################################################################### # Subsurface Model # ---------------- +loop_current = 1.0 inp = { - 'res': np.array([2e14, 1.0, 100]), - 'times': times, - 'inp': { - 'src': src_dipole, - 'rec': rec, - 'depth': [0.0, 50.0], - 'verb': 1, - } + "src": src_dipole, + "rec": rec, + "depth": [0.0, 50.0], + "freqtime": times, + "signal": {"nodes": nodes, "amplitudes": amplitudes, "signal": -1}, + "srcpts": 10, + "msrc": False, + "ftarg": {"dlf": "key_81_2009"}, + "htarg": {"dlf": "key_101_2009", "pts_per_dec": -1}, + "verb": 1, } -apply_cole_cole = { - 'm': np.array([0.0, 0.05, 0.15]), - 'tau': np.array([0.0, 0.005, 0.02]), - 'c': np.array([0.0, 0.4, 0.6]), +param_ip = { + "m": np.array([0.0, 0.05, 0.15]), + "tau": np.array([0.0, 0.005, 0.02]), + "c": np.array([0.0, 0.4, 0.6]), } -apply_vrm = { - 'dchi': np.array([0.0, 0.005, 0.02]), - 'mu': np.array([mu_0, mu_0 * 1.05, mu_0]), +param_vrm = { + "dchi": np.array([0.0, 0.005, 0.02]), + "mu": np.array([mu_0, mu_0 * 1.05, mu_0]), } +# Build empymod model dict with VRM and/or Cole-Cole hooks as needed +res = np.array([2e14, 1.0, 100]) +res_vrm = {"res": res, "func_zeta": vrm_from_mu, **param_vrm} +res_ip = {"res": res, "func_eta": pelton_cole_cole_model, **param_ip} +res_both = {**res_vrm, **res_ip} + ############################################################################### # Compute responses @@ -445,19 +239,23 @@ def simulate_empymod(inp, res, times, compute_B_field=True, loop_current=1.0, for compute_B_field in [True, False]: - inp['compute_B_field'] = compute_B_field + # Source strength and receiver type depends on output field type: + # - For B field: `mrec=True` yields H field; mu_0 * current yields B=μH in + # time domain. + # - For dB/dt: `mrec='b'` multiplies by iωμ in frequency domain. + inp["mrec"] = True if compute_B_field else "b" + inp["strength"] = mu_0 * loop_current if compute_B_field else loop_current out = {} # Conductivity - out['σ'] = simulate_empymod(**inp) + out["σ"] = empymod.bipole(**inp, res=res).sum(axis=-1) # IP - out['IP'] = simulate_empymod(**inp, apply_cole_cole=apply_cole_cole,) + out["IP"] = empymod.bipole(**inp, res=res_ip).sum(axis=-1) # VRM - out['VRM'] = simulate_empymod(**inp, apply_vrm=apply_vrm) + out["VRM"] = empymod.bipole(**inp, res=res_vrm).sum(axis=-1) # IP + VRM - out['IP+VRM'] = simulate_empymod( - **inp, apply_cole_cole=apply_cole_cole, apply_vrm=apply_vrm) + out["IP+VRM"] = empymod.bipole(**inp, res=res_both).sum(axis=-1) results["B field" if compute_B_field else "dB/dt"] = out @@ -467,10 +265,10 @@ def simulate_empymod(inp, res, times, compute_B_field=True, loop_current=1.0, # ------------ fig, axs = plt.subplots( - 2, 2, figsize=(10, 7), sharex=True, sharey=True, layout='constrained' + 2, 2, figsize=(10, 7), sharex=True, sharey="row", layout="constrained" ) -in_out = ['Outside', 'Inside'] +in_out = ["Outside", "Inside"] # Loop over B-field - dB/dt for i, compute_B_field in enumerate(["B field", "dB/dt"]): @@ -481,17 +279,17 @@ def simulate_empymod(inp, res, times, compute_B_field=True, loop_current=1.0, # Loop over cases for k, v in results[compute_B_field].items(): - axs[i, ii].loglog(times*1e3, np.abs(v[ii, :]), label=k) + axs[i, ii].loglog(times*1e3, np.abs(v[:, ii]), label=k) axs[i, ii].set_title(f"{compute_B_field}, {in_out[ii]}") axs[i, ii].legend() for ax in axs[1, :]: - ax.set_xlabel('Time (ms)') + ax.set_xlabel("Time (ms)") axs[0, 0].set_ylabel("$B_z$ (T)") axs[1, 0].set_ylabel("$dB_z/dt$ (T/s)") -plt.suptitle("Comparison of σ, IP, VRM, and IP+VRM", fontsize=14) +plt.suptitle("Comparison of σ, IP, VRM, and IP+VRM", fontsize=18) plt.show() ############################################################################### diff --git a/examples/time_domain/tem_temfast.py b/examples/time_domain/tem_temfast.py index c221ef2..2744ea8 100644 --- a/examples/time_domain/tem_temfast.py +++ b/examples/time_domain/tem_temfast.py @@ -5,10 +5,10 @@ **In this example we compute the TEM response from the TEM-FAST 48 system.** This example was contributed by Lukas Aigner (`@aignerlukas -`_), who was interested -in modelling the TEM-FAST system, which is used at the TU Wien. -If you are interested and want to use this work please have a look at the -corresponding paper Aigner et al. (2024). +`_), who was interested in modelling the +TEM-FAST system, which is used at the TU Wien. If you are interested and want +to use this work please have a look at the corresponding paper Aigner et al. +(2024). The modeller ``empymod`` models the electromagnetic (EM) full wavefield Greens function for electric and magnetic point sources and receivers. As such, it can @@ -16,12 +16,8 @@ particular EM method and survey layout can be tricky, as there are many more things involved than just computing the EM Greens function. -What is not included in ``empymod`` at this moment (but hopefully in the -future), but is required to model TEM data, is to **account for arbitrary -source waveform**, and to apply a **lowpass filter**. So we generate these two -things here, and create our own wrapper to model TEM data. See also the example -:ref:`sphx_glr_gallery_tdomain_tem_walktem.py`, on which this example -builds upon. +See also the example :ref:`sphx_glr_gallery_tdomain_tem_walktem.py`, on which +this example builds upon. **References** @@ -36,10 +32,9 @@ import empymod import numpy as np import matplotlib.pyplot as plt -from scipy.special import roots_legendre -from matplotlib.ticker import LogLocator, NullFormatter -from scipy.interpolate import InterpolatedUnivariateSpline as iuSpline + plt.style.use('ggplot') + # sphinx_gallery_thumbnail_number = 2 ############################################################################### @@ -49,257 +44,99 @@ # The TEM-FASt system uses a "time-key" value to determine the number of gates, # the front ramp and the length of the current pulse. # We are using values that correspond to a time-key of 5. -turn_on_ramp = -3.0E-06 -turn_off_ramp = 0.95E-06 -on_time = 3.75E-03 +turn_on_ramp = -3e-6 +turn_off_ramp = 0.95e-6 +on_time = 3.75e-3 injected_current = 4.1 # Ampere -time_gates = np.r_[4.060e+00, 5.070e+00, 6.070e+00, 7.080e+00, - 8.520e+00, 1.053e+01, 1.255e+01, 1.456e+01, - 1.744e+01, 2.146e+01, 2.549e+01, 2.950e+01, - 3.528e+01, 4.330e+01, 5.140e+01, 5.941e+01, # time-key 1 - 7.160e+01, 8.760e+01, 1.036e+02, 1.196e+02, # time-key 2 - 1.436e+02, 1.756e+02, 2.076e+02, 2.396e+02, # time-key 3 - 2.850e+02, 3.500e+02, 4.140e+02, 4.780e+02, # time-key 4 - 5.700e+02, 6.990e+02, 8.280e+02, 9.560e+02, # time-key 5 - ] * 1e-6 # from us to s - -waveform_times = np.r_[turn_on_ramp - on_time, -on_time, - 0.000E+00, turn_off_ramp] -waveform_current = np.r_[0.0, injected_current, injected_current, 0.0] - -plt.figure() -plt.title('Waveform') -plt.plot(np.r_[-9, waveform_times*1e3, 2], np.r_[0, waveform_current, 0]) -plt.xlabel('Time (ms)') -plt.xlim([-4, 0.5]) +time_gates = np.array([ + 4.060e+00, 5.070e+00, 6.070e+00, 7.080e+00, + 8.520e+00, 1.053e+01, 1.255e+01, 1.456e+01, + 1.744e+01, 2.146e+01, 2.549e+01, 2.950e+01, + 3.528e+01, 4.330e+01, 5.140e+01, 5.941e+01, # time-key 1 + 7.160e+01, 8.760e+01, 1.036e+02, 1.196e+02, # time-key 2 + 1.436e+02, 1.756e+02, 2.076e+02, 2.396e+02, # time-key 3 + 2.850e+02, 3.500e+02, 4.140e+02, 4.780e+02, # time-key 4 + 5.700e+02, 6.990e+02, 8.280e+02, 9.560e+02, # time-key 5 +]) * 1e-6 # from us to s + +nodes = np.array([turn_on_ramp - on_time, -on_time, 0, turn_off_ramp]) +amplitudes = np.array([0.0, injected_current, injected_current, 0.0]) + +fig, ax = plt.subplots(1, 1, constrained_layout=True) +ax.set_title('Waveform') +ax.plot(np.r_[-5, nodes*1e3, 1.5], np.r_[0, amplitudes, 0]) +ax.set_xlabel('Time (ms)') +ax.set_xlim([-5, 1.5]) ############################################################################### # 2. ``empymod`` implementation # ----------------------------- -def waveform(times, resp, times_wanted, wave_time, wave_amp, nquad=3): - """Apply a source waveform to the signal. - - Parameters - ---------- - times : ndarray - Times of computed input response; should start before and end after - `times_wanted`. - - resp : ndarray - EM-response corresponding to `times`. - - times_wanted : ndarray - Wanted times. - - wave_time : ndarray - Time steps of the wave. - - wave_amp : ndarray - Amplitudes of the wave corresponding to `wave_time`, usually - in the range of [0, 1]. - - nquad : int - Number of Gauss-Legendre points for the integration. Default is 3. - - Returns - ------- - resp_wanted : ndarray - EM field for `times_wanted`. - - """ - - # Interpolate on log. - PP = iuSpline(np.log10(times), resp) - - # Wave time steps. - dt = np.diff(wave_time) - dI = np.diff(wave_amp) - dIdt = dI/dt - - # Gauss-Legendre Quadrature; 3 is generally good enough. - # (Roots/weights could be cached.) - g_x, g_w = roots_legendre(nquad) - - # Pre-allocate output. - resp_wanted = np.zeros_like(times_wanted) - - # Loop over wave segments. - for i, cdIdt in enumerate(dIdt): - - # We only have to consider segments with a change of current. - if cdIdt == 0.0: - continue - - # If wanted time is before a wave element, ignore it. - ind_a = wave_time[i] < times_wanted - if ind_a.sum() == 0: - continue - - # If wanted time is within a wave element, we cut the element. - ind_b = wave_time[i+1] > times_wanted[ind_a] - - # Start and end for this wave-segment for all times. - ta = times_wanted[ind_a]-wave_time[i] - tb = times_wanted[ind_a]-wave_time[i+1] - tb[ind_b] = 0.0 # Cut elements - - # Gauss-Legendre for this wave segment. See - # https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval - # for the change of interval, which makes this a bit more complex. - logt = np.log10(np.outer((tb-ta)/2, g_x)+(ta+tb)[:, None]/2) - fact = (tb-ta)/2*cdIdt - resp_wanted[ind_a] += fact*np.sum(np.array(PP(logt)*g_w), axis=1) - - return resp_wanted - - -############################################################################### -def get_time(time, r_time): - """Additional time for ramp. - - Because of the arbitrary waveform, we need to compute some times before and - after the actually wanted times for interpolation of the waveform. - - Some implementation details: The actual times here don't really matter. We - create a vector of time.size+2, so it is similar to the input times and - accounts that it will require a bit earlier and a bit later times. Really - important are only the minimum and maximum times. The Fourier DLF, with - `pts_per_dec=-1`, computes times from minimum to at least the maximum, - where the actual spacing is defined by the filter spacing. It subsequently - interpolates to the wanted times. Afterwards, we interpolate those again to - compute the actual waveform response. - - Note: We could first call `waveform`, and get the actually required times - from there. This would make this function obsolete. It would also - avoid the double interpolation, first in `empymod.model.time` for the - Fourier DLF with `pts_per_dec=-1`, and second in `waveform`. Doable. - Probably not or marginally faster. And the code would become much - less readable. - - Parameters - ---------- - time : ndarray - Desired times - - r_time : ndarray - Waveform times +# +# Here we collect the necessary input for empymod to model temfast. - Returns - ------- - time_req : ndarray - Required times - """ - tmin = np.log10(max(time.min()-r_time.max(), 1e-10)) - tmax = np.log10(time.max()-r_time.min()) - return np.logspace(tmin, tmax, time.size+2) +def bandpass(inp, p_dict): + """Butterworth-type filter (implemented from simpegEM1D.Waveforms.py).""" + cutofffreq = 1e8 # Determined empirically for TEM-FAST + h = (1 + 1j*p_dict["freq"]/cutofffreq)**-1 + h *= (1 + 1j*p_dict["freq"]/3e5)**-1 + p_dict["EM"] *= h[:, None] -############################################################################### -def temfast(off_time, waveform_times, model, square_side=12.5): +def temfast(depth, res): """Custom wrapper of empymod.model.bipole. Here, we compute TEM-FAST data using the ``empymod.model.bipole`` routine - as an example. This function is based upon the Walk TEM example. - - We model the big source square loop by computing only half of one side of - the electric square loop and approximating the finite length dipole with 3 - point dipole sources. The result is then multiplied by 8, to account for - all eight half-sides of the square loop. - - The implementation here assumes a central loop configuration, where the - receiver (1 m2 area) is at the origin, and the source is a - square_side x square_side m electric loop, centered around the origin. - - Note: This approximation of only using half of one of the four sides - obviously only works for central, horizontal square loops. If your - loop is arbitrary rotated, then you have to model all four sides of - the loop and sum it up. + as an example. Everything is fixed except for the depth and the resistivity + models. Parameters ---------- - off_time : ndarray - times at which the secondary magnetic field will be measured - - waveform_times : ndarray - Depths of the resistivity model (see ``empymod.model.bipole`` for more - info.) - depth : ndarray - Depths of the resistivity model (see ``empymod.model.bipole`` for more - info.) + Depths of the resistivity model interfaces (see + ``empymod.model.bipole`` for more info). res : ndarray Resistivities of the resistivity model (see ``empymod.model.bipole`` - for more info.) - - square_side : float - sige length of the square loop in meter. + for more info). Returns ------- - TEM-FAST waveform : EMArray - TEM-FAST response (dB/dt). + TEM-FAST : EMArray + TEM-FAST response [dB/dt]. """ - if 'm' in model: - depth = model['depth'] - res = model - del res['depth'] - else: - res = model['res'] - depth = model['depth'] - - # === GET REQUIRED TIMES === - time = get_time(off_time, waveform_times) - - # === GET REQUIRED FREQUENCIES === - time, freq, ft, ftarg = empymod.utils.check_time( - time=time, # Required times - signal=1, # Switch-on response - ft='dlf', # Use DLF - ftarg={'dlf': 'key_601_2009'}, - verb=2, - ) + # The waveform signal + signal = {'nodes': nodes, 'amplitudes': amplitudes, 'signal': 1} - # === COMPUTE FREQUENCY-DOMAIN RESPONSE === + # === COMPUTE RESPONSE === # We only define a few parameters here. You could extend this for any # parameter possible to provide to empymod.model.bipole. + + # We model the square with 1/2 of one side. This makes it faster, but it + # will only work for a horizontal square loop, with a central receiver. + square_side = 12.5 hs = square_side / 2 # half side length + EM = empymod.model.bipole( - src=[hs, hs, 0, hs, 0, 0], # El. dipole source; half of one side. - rec=[0, 0, 0, 0, 90], # Receiver at the origin, vertical. - depth=depth, # Depth-model, including air-interface. - res=res, # if with IP, res is a dictionary with - # all params and the function - freqtime=freq, # Required frequencies. - mrec=True, # It is an el. source, but a magn. rec. - strength=8, # To account for 4 sides of square loop. - srcpts=3, # Approx. the finite dip. with 3 points. - htarg={'dlf': 'key_401_2009'}, # Short filter, so fast. + src=[hs, hs, 0, hs, 0, 0], # El. dipole source; half of one side. + rec=[0, 0, 0, 0, 90], # Receiver at the origin, vertical. + depth=depth, # Depth-model. + res=res, # Resistivity model. + freqtime=time_gates, # Wanted times. + signal=signal, # Waveform + mrec="b", # Receiver: dB/dt + strength=8, # To account for 8 quarters of square. + srcpts=3, # Approx. the finite dip. with 3 points. + ftarg={"dlf": "key_81_2009"}, # Shorter, faster filters. + htarg={"dlf": "key_101_2009", "pts_per_dec": -1}, + bandpass={"func": bandpass}, ) - # Multiply the frequecny-domain result with - # \mu for H->B, and i\omega for B->dB/dt. - EM *= 2j*np.pi*freq*4e-7*np.pi - - # === Butterworth-type filter (implemented from simpegEM1D.Waveforms.py)=== - cutofffreq = 1e8 # determined empirically for TEM-FAST - h = (1+1j*freq/cutofffreq)**-1 # First order type - h *= (1+1j*freq/3e5)**-1 - EM *= h - - # === CONVERT TO TIME DOMAIN === - delay_rst = 0 # unknown for TEM-FAST, therefore 0 - EM, _ = empymod.model.tem(EM[:, None], np.array([1]), - freq, time+delay_rst, 1, ft, ftarg) - EM = np.squeeze(EM) - - # === APPLY WAVEFORM === - return waveform(time, EM, off_time, waveform_times, waveform_current) + return EM ############################################################################### @@ -313,7 +150,7 @@ def pelton_res(inp, p_dict): # print('\n shape: p_dict["freq"]\n', p_dict['freq'].shape) iwtc = np.outer(2j*np.pi*p_dict['freq'], inp['tau'])**inp['c'] - rhoH = inp['rho_0'] * (1 - inp['m']*(1 - 1/(1 + iwtc))) + rhoH = inp['res'] * (1 - inp['m']*(1 - 1/(1 + iwtc))) rhoV = rhoH*p_dict['aniso']**2 # Add electric permittivity contribution @@ -324,77 +161,50 @@ def pelton_res(inp, p_dict): ############################################################################### -# 3. Computation non-IP -# --------------------- - -depths = [8, 20] -rhos = [25, 5, 50] -model = {'depth': np.r_[0, depths], - 'res': np.r_[2e14, rhos]} +# 3. Computate responses +# ---------------------- -# Compute conductive model -response = temfast(off_time=time_gates, waveform_times=waveform_times, - model=model) +depth = [0, 8, 20] +rhos = [2e14, 25, 5, 50] +rhos_ip = { + 'res': rhos, + 'm': np.array([0, 0, 0.9, 0]), + 'tau': np.array([1e-7, 1e-6, 5e-4, 1e-6]), + 'c': np.array([0.01, 0, 0.9, 0]), + 'func_eta': pelton_res, +} -############################################################################### -# 4. Computation with IP -# ---------------------- -depths = [8, 20] -rhos = [25, 5, 50] -charg = np.r_[0, 0.9, 0] -taus = np.r_[1e-6, 5e-4, 1e-6] -cs = np.r_[0, 0.9, 0] - -eta_func = pelton_res -depth = np.r_[0, depths] -model = {'depth': depth, - 'res': np.r_[2e14, rhos], - 'rho_0': np.r_[2e14, rhos], - 'm': np.r_[0, charg], - 'tau': np.r_[1e-7, taus], - 'c': np.r_[0.01, cs], - 'func_eta': eta_func} - +# Compute conductive model +response = temfast(depth=depth, res=rhos) # Compute conductive model -response_ip = temfast(off_time=time_gates, waveform_times=waveform_times, - model=model) +response_ip = temfast(depth=depth, res=rhos_ip) ############################################################################### -# 5. Comparison +# 4. Comparison # ------------- -plt.figure(figsize=(5, 5), constrained_layout=True) +fig, ax = plt.subplots(1, 1, constrained_layout=True) -# Plot result of model 1 -ax1 = plt.subplot(111) -plt.title('TEM-FAST responses') +ax.set_title('TEM-FAST responses') # empymod -plt.plot(time_gates, response, 'r.--', ms=7, label="response") -plt.plot(time_gates, abs(response_ip), 'kx:', ms=7, label="response with IP") +ax.loglog(time_gates, response, 'r.--', ms=7, label="response") +ax.loglog(time_gates, abs(response_ip), 'kx:', ms=7, label="response with IP") sub0 = response_ip[response_ip < 0] tg_sub0 = time_gates[response_ip < 0] -plt.plot(tg_sub0, abs(sub0), marker='s', ls='none', - mfc='none', ms=8, mew=1, - mec='c', label="negative readings") +ax.loglog(tg_sub0, abs(sub0), marker='s', ls='none', mfc='none', + ms=8, mew=1, mec='c', label="negative readings") # Plot settings -plt.xscale('log') -plt.yscale('log') -plt.xlabel("Time(s)") -plt.ylabel(r"$\mathrm{d}\mathrm{B}_\mathrm{z}\,/\,\mathrm{d}t$") -plt.grid(which='both', c='w') -plt.legend(title='Data', loc=1) - - -# Force minor ticks on logscale -ax1.yaxis.set_minor_locator(LogLocator(subs='all', numticks=20)) -ax1.yaxis.set_minor_formatter(NullFormatter()) -plt.grid(which='both', c='w') +ax.set_xlabel("Time(s)") +ax.set_ylabel(r"$\mathrm{d}\mathrm{B}_\mathrm{z}\,/\,\mathrm{d}t$") +ax.grid(which='both') +ax.legend() + ############################################################################### empymod.Report() diff --git a/examples/time_domain/tem_walktem.py b/examples/time_domain/tem_walktem.py index a6e09ac..b0d0abd 100644 --- a/examples/time_domain/tem_walktem.py +++ b/examples/time_domain/tem_walktem.py @@ -11,13 +11,8 @@ **In this example we are going to compute a TEM response, in particular from the system** WalkTEM, and compare it with data obtained from `AarhusInv `_. However, you can use and adapt this -example to model other TEM systems, such as skyTEM, SIROTEM, TEM-FAST, or any -other system. - -What is not included in ``empymod`` at this moment (but hopefully in the -future), but is required to model TEM data, is to **account for arbitrary -source waveform**, and to apply a **lowpass filter**. So we generate these two -things here, and create our own wrapper to model TEM data. +example to model other TEM systems, such as skyTEM, SIROTEM, TEM-FAST +(:ref:`sphx_glr_gallery_tdomain_tem_temfast.py`), or any other system. The incentive for this example came from Leon Foks (`@leonfoks `_) for `GeoBIPy @@ -29,10 +24,9 @@ import empymod import numpy as np import matplotlib.pyplot as plt -from scipy.special import roots_legendre -from matplotlib.ticker import LogLocator, NullFormatter -from scipy.interpolate import InterpolatedUnivariateSpline as iuSpline + plt.style.use('ggplot') + # sphinx_gallery_thumbnail_number = 2 ############################################################################### @@ -107,267 +101,132 @@ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Low moment -lm_waveform_times = np.r_[-1.041E-03, -9.850E-04, 0.000E+00, 4.000E-06] -lm_waveform_current = np.r_[0.0, 1.0, 1.0, 0.0] +lm_waveform_times = np.array([-1.041e-3, -9.850e-4, 0, 4e-6]) +lm_waveform_current = np.array([0.0, 1.0, 1.0, 0.0]) # High moment -hm_waveform_times = np.r_[-8.333E-03, -8.033E-03, 0.000E+00, 5.600E-06] -hm_waveform_current = np.r_[0.0, 1.0, 1.0, 0.0] +hm_waveform_times = np.array([-8.333e-3, -8.033e-3, 0, 5.6e-6]) +hm_waveform_current = np.array([0.0, 1.0, 1.0, 0.0]) -plt.figure() -plt.title('Waveforms') -plt.plot(np.r_[-9, lm_waveform_times*1e3, 2], np.r_[0, lm_waveform_current, 0], - label='Low moment') -plt.plot(np.r_[-9, hm_waveform_times*1e3, 2], np.r_[0, hm_waveform_current, 0], - '-.', label='High moment') -plt.xlabel('Time (ms)') -plt.xlim([-9, 0.5]) -plt.legend() +# Plot them +fig, ax = plt.subplots(1, 1, constrained_layout=True) +ax.set_title('Waveforms') +ax.plot(np.r_[-9, lm_waveform_times*1e3, 2], np.r_[0, lm_waveform_current, 0], + label='Low moment') +ax.plot(np.r_[-9, hm_waveform_times*1e3, 2], np.r_[0, hm_waveform_current, 0], + '-.', label='High moment') +ax.set_xlabel('Time (ms)') +ax.set_xlim([-9, 0.5]) +ax.legend() ############################################################################### # 2. ``empymod`` implementation # ----------------------------- -def waveform(times, resp, times_wanted, wave_time, wave_amp, nquad=3): - """Apply a source waveform to the signal. - - Parameters - ---------- - times : ndarray - Times of computed input response; should start before and end after - `times_wanted`. - - resp : ndarray - EM-response corresponding to `times`. - - times_wanted : ndarray - Wanted times. - - wave_time : ndarray - Time steps of the wave. - - wave_amp : ndarray - Amplitudes of the wave corresponding to `wave_time`, usually - in the range of [0, 1]. - - nquad : int - Number of Gauss-Legendre points for the integration. Default is 3. - - Returns - ------- - resp_wanted : ndarray - EM field for `times_wanted`. - - """ - - # Interpolate on log. - PP = iuSpline(np.log10(times), resp) - - # Wave time steps. - dt = np.diff(wave_time) - dI = np.diff(wave_amp) - dIdt = dI/dt - - # Gauss-Legendre Quadrature; 3 is generally good enough. - # (Roots/weights could be cached.) - g_x, g_w = roots_legendre(nquad) - - # Pre-allocate output. - resp_wanted = np.zeros_like(times_wanted) - - # Loop over wave segments. - for i, cdIdt in enumerate(dIdt): - - # We only have to consider segments with a change of current. - if cdIdt == 0.0: - continue - - # If wanted time is before a wave element, ignore it. - ind_a = wave_time[i] < times_wanted - if ind_a.sum() == 0: - continue - - # If wanted time is within a wave element, we cut the element. - ind_b = wave_time[i+1] > times_wanted[ind_a] - - # Start and end for this wave-segment for all times. - ta = times_wanted[ind_a]-wave_time[i] - tb = times_wanted[ind_a]-wave_time[i+1] - tb[ind_b] = 0.0 # Cut elements - - # Gauss-Legendre for this wave segment. See - # https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval - # for the change of interval, which makes this a bit more complex. - logt = np.log10(np.outer((tb-ta)/2, g_x)+(ta+tb)[:, None]/2) - fact = (tb-ta)/2*cdIdt - resp_wanted[ind_a] += fact*np.sum(np.array(PP(logt)*g_w), axis=1) - - return resp_wanted - - -############################################################################### -def get_time(time, r_time): - """Additional time for ramp. - - Because of the arbitrary waveform, we need to compute some times before and - after the actually wanted times for interpolation of the waveform. - - Some implementation details: The actual times here don't really matter. We - create a vector of time.size+2, so it is similar to the input times and - accounts that it will require a bit earlier and a bit later times. Really - important are only the minimum and maximum times. The Fourier DLF, with - `pts_per_dec=-1`, computes times from minimum to at least the maximum, - where the actual spacing is defined by the filter spacing. It subsequently - interpolates to the wanted times. Afterwards, we interpolate those again to - compute the actual waveform response. - - Note: We could first call `waveform`, and get the actually required times - from there. This would make this function obsolete. It would also - avoid the double interpolation, first in `empymod.model.time` for the - Fourier DLF with `pts_per_dec=-1`, and second in `waveform`. Doable. - Probably not or marginally faster. And the code would become much - less readable. - - Parameters - ---------- - time : ndarray - Desired times - - r_time : ndarray - Waveform times - - Returns - ------- - time_req : ndarray - Required times - """ - tmin = np.log10(max(time.min()-r_time.max(), 1e-10)) - tmax = np.log10(time.max()-r_time.min()) - return np.logspace(tmin, tmax, time.size+2) +# +# We model the big source square loop by computing only half of one side of +# the electric square loop and approximating the finite length dipole with 3 +# point dipole sources. The result is then multiplied by 8, to account for +# all eight half-sides of the square loop. +# +# The implementation here assumes a central loop configuration, where the +# receiver (1 m² area) is at the origin, and the source is a 40x40 m electric +# loop, centered around the origin. +# +# Note: This approximation of only using half of one of the four sides +# obviously only works for central, horizontal square loops. If your loop +# is arbitrary rotated, then you have to model all four sides of the loop +# and sum it up. +# +# As an example, if the receiver wouldn't be in the center, we would have to +# model the actual complete loop (no symmetry to take advantage of). +# +# .. code-block:: python +# +# EM = empymod.model.bipole( +# src=[[20, 20, -20, -20], # x1 +# [20, -20, -20, 20], # x2 +# [-20, 20, 20, -20], # y1 +# [20, 20, -20, -20], # y2 +# 0, 0], # z1, z2 +# strength=1, +# # ... all other parameters remain the same +# ) +# EM = EM.sum(axis=1) # Sum all source dipoles + +def bandpass(inp, p_dict): + """Butterworth-type filter (implemented from simpegEM1D.Waveforms.py).""" + cutofffreq = 4.5e5 # As stated in the WalkTEM manual + h = (1 + 1j*p_dict["freq"]/cutofffreq)**-1 + h *= (1 + 1j*p_dict["freq"]/3e5)**-1 + p_dict["EM"] *= h[:, None] -############################################################################### def walktem(moment, depth, res): """Custom wrapper of empymod.model.bipole. Here, we compute WalkTEM data using the ``empymod.model.bipole`` routine as - an example. We could achieve the same using ``empymod.model.dipole`` or - ``empymod.model.loop``. - - We model the big source square loop by computing only half of one side of - the electric square loop and approximating the finite length dipole with 3 - point dipole sources. The result is then multiplied by 8, to account for - all eight half-sides of the square loop. - - The implementation here assumes a central loop configuration, where the - receiver (1 m2 area) is at the origin, and the source is a 40x40 m electric - loop, centered around the origin. - - Note: This approximation of only using half of one of the four sides - obviously only works for central, horizontal square loops. If your - loop is arbitrary rotated, then you have to model all four sides of - the loop and sum it up. + an example. Everything is fixed except for the moment, the depth model, and + the resistivity model. Parameters ---------- - moment : str {'lm', 'hm'} - Moment. If 'lm', above defined ``lm_off_time``, ``lm_waveform_times``, + moment : str {"lm", "hm"} + Moment. If "lm", above defined ``lm_off_time``, ``lm_waveform_times``, and ``lm_waveform_current`` are used. Else, the corresponding ``hm_``-parameters. depth : ndarray - Depths of the resistivity model (see ``empymod.model.bipole`` for more - info.) + Depths of the resistivity model interfaces (see + ``empymod.model.bipole`` for more info), without 0. res : ndarray Resistivities of the resistivity model (see ``empymod.model.bipole`` - for more info.) + for more info), without air. Returns ------- WalkTEM : EMArray - WalkTEM response (dB/dt). + WalkTEM response [dB/dt]. """ - # Get the measurement time and the waveform corresponding to the provided - # moment. - if moment == 'lm': + # Get measurement time and waveform corresponding to the provided moment. + if moment == "lm": off_time = lm_off_time - waveform_times = lm_waveform_times - waveform_current = lm_waveform_current - elif moment == 'hm': + nodes = lm_waveform_times + amplitudes = lm_waveform_current + elif moment == "hm": off_time = hm_off_time - waveform_times = hm_waveform_times - waveform_current = hm_waveform_current + nodes = hm_waveform_times + amplitudes = hm_waveform_current else: raise ValueError("Moment must be either 'lm' or 'hm'!") - # === GET REQUIRED TIMES === - time = get_time(off_time, waveform_times) + # Collect signal + signal = {"nodes": nodes, "amplitudes": amplitudes, "signal": 1} + delay = 1.8e-7 # As stated in the WalkTEM manual - # === GET REQUIRED FREQUENCIES === - time, freq, ft, ftarg = empymod.utils.check_time( - time=time, # Required times - signal=1, # Switch-on response - ft='dlf', # Use DLF - ftarg={'dlf': 'key_81_2009'}, # Short, fast filter; if you - verb=2, # need higher accuracy choose a longer filter. - ) - - # === COMPUTE FREQUENCY-DOMAIN RESPONSE === + # === COMPUTE RESPONSE === # We only define a few parameters here. You could extend this for any # parameter possible to provide to empymod.model.bipole. EM = empymod.model.bipole( - src=[20, 20, 0, 20, 0, 0], # El. dipole source; half of one side. - rec=[0, 0, 0, 0, 90], # Receiver at the origin, vertical. - depth=np.r_[0, depth], # Depth-model, adding air-interface. - res=np.r_[2e14, res], # Provided resistivity model, adding air. - # aniso=aniso, # Here you could implement anisotropy... - # # ...or any parameter accepted by bipole. - freqtime=freq, # Required frequencies. - mrec=True, # It is an el. source, but a magn. rec. - strength=8, # To account for 4 sides of square loop. - srcpts=3, # Approx. the finite dip. with 3 points. - htarg={'dlf': 'key_101_2009'}, # Short filter, so fast. + src=[20, 20, 0, 20, 0, 0], # El. dipole source; half of one side. + rec=[0, 0, 0, 0, 90], # Receiver at the origin, vertical. + depth=np.r_[0, depth], # Depth-model, adding air-interface. + res=np.r_[2e14, res], # Provided resistivity model, adding air. + freqtime=off_time + delay, # Wanted times + signal=signal, # Waveform + mrec="b", # Receiver: dB/dt + strength=8, # To account for 8 quarters of square. + srcpts=3, # Approx. the finite dip. with 3 points. + ftarg={"dlf": "key_81_2009"}, # Shorter, faster filters. + htarg={"dlf": "key_101_2009", "pts_per_dec": -1}, + bandpass={"func": bandpass}, ) - # Note: If the receiver wouldn't be in the center, we would have to model - # the actual complete loop (no symmetry to take advantage of). - # - # EM = empymod.model.bipole( - # src=[[20, 20, -20, -20], # x1 - # [20, -20, -20, 20], # x2 - # [-20, 20, 20, -20], # y1 - # [20, 20, -20, -20], # y2 - # 0, 0], # z1, z2 - # strength=1, - # # ... all other parameters remain the same - # ) - # EM = EM.sum(axis=1) # Sum all source dipoles - - # Multiply the frequecny-domain result with - # \mu for H->B, and i\omega for B->dB/dt. - EM *= 2j*np.pi*freq*4e-7*np.pi - - # === Butterworth-type filter (implemented from simpegEM1D.Waveforms.py)=== - # Note: Here we just apply one filter. But it seems that WalkTEM can apply - # two filters, one before and one after the so-called front gate - # (which might be related to ``delay_rst``, I am not sure about that - # part.) - cutofffreq = 4.5e5 # As stated in the WalkTEM manual - h = (1+1j*freq/cutofffreq)**-1 # First order type - h *= (1+1j*freq/3e5)**-1 - EM *= h - - # === CONVERT TO TIME DOMAIN === - delay_rst = 1.8e-7 # As stated in the WalkTEM manual - EM, _ = empymod.model.tem(EM[:, None], np.array([1]), - freq, time+delay_rst, 1, ft, ftarg) - EM = np.squeeze(EM) - - # === APPLY WAVEFORM === - return waveform(time, EM, off_time, waveform_times, waveform_current) + + return EM ############################################################################### @@ -382,71 +241,61 @@ def walktem(moment, depth, res): lm_empymod_con = walktem('lm', depth=[30], res=[10, 1]) hm_empymod_con = walktem('hm', depth=[30], res=[10, 1]) + ############################################################################### # 4. Comparison # ------------- -plt.figure(figsize=(9, 5), constrained_layout=True) +fig, axs = plt.subplots(1, 2, figsize=(9, 5), constrained_layout=True) +ax1, ax2 = axs # Plot result resistive model -ax1 = plt.subplot(121) -plt.title('Resistive Model') +ax1.set_title("Resistive Model") # AarhusInv -plt.plot(lm_off_time, lm_aarhus_res, 'd', mfc='.4', mec='.4', - label="Aarhus LM") -plt.plot(hm_off_time, hm_aarhus_res, 's', mfc='.4', mec='.4', - label="Aarhus HM") +ax1.loglog(lm_off_time, lm_aarhus_res, "d", c=".4", label="Aarhus LM") +ax1.loglog(hm_off_time, hm_aarhus_res, "s", c=".4", label="Aarhus HM") # empymod -plt.plot(lm_off_time, lm_empymod_res, 'r+', ms=7, label="empymod LM") -plt.plot(hm_off_time, hm_empymod_res, 'cx', label="empymod HM") +ax1.loglog(lm_off_time, lm_empymod_res, "r+", ms=7, label="empymod LM") +ax1.loglog(hm_off_time, hm_empymod_res, "cx", label="empymod HM") # Difference -plt.plot(lm_off_time, np.abs((lm_aarhus_res - lm_empymod_res)), 'm.') -plt.plot(hm_off_time, np.abs((hm_aarhus_res - hm_empymod_res)), 'b.') +ax1.loglog(lm_off_time, np.abs((lm_aarhus_res - lm_empymod_res)), "m.") +ax1.loglog(hm_off_time, np.abs((hm_aarhus_res - hm_empymod_res)), "b.") -# Plot settings -plt.xscale('log') -plt.yscale('log') -plt.xlabel("Time(s)") -plt.ylabel(r"$\mathrm{d}\mathrm{B}_\mathrm{z}\,/\,\mathrm{d}t$") -plt.grid(which='both', c='w') -plt.legend(title='Data', loc=1) +# Legend +ax1.legend(title="Data") # Plot result conductive model -ax2 = plt.subplot(122) -plt.title('Conductive Model') +ax2.set_title("Conductive Model") ax2.yaxis.set_label_position("right") ax2.yaxis.tick_right() # AarhusInv -plt.plot(lm_off_time, lm_aarhus_con, 'd', mfc='.4', mec='.4') -plt.plot(hm_off_time, hm_aarhus_con, 's', mfc='.4', mec='.4') +ax2.loglog(lm_off_time, lm_aarhus_con, "d", c=".4") +ax2.loglog(hm_off_time, hm_aarhus_con, "s", c=".4") # empymod -plt.plot(lm_off_time, lm_empymod_con, 'r+', ms=7) -plt.plot(hm_off_time, hm_empymod_con, 'cx') +ax2.loglog(lm_off_time, lm_empymod_con, "r+", ms=7) +ax2.loglog(hm_off_time, hm_empymod_con, "cx") # Difference -plt.plot(lm_off_time, np.abs((lm_aarhus_con - lm_empymod_con)), 'm.', - label=r"$|\Delta_\mathrm{LM}|$") -plt.plot(hm_off_time, np.abs((hm_aarhus_con - hm_empymod_con)), 'b.', - label=r"$|\Delta_\mathrm{HM}|$") - -# Plot settings -plt.xscale('log') -plt.yscale('log') -plt.xlabel("Time(s)") -plt.ylabel(r"$\mathrm{d}\mathrm{B}_\mathrm{z}\,/\,\mathrm{d}t$") -plt.legend(title='Difference', loc=3) - -# Force minor ticks on logscale -ax1.yaxis.set_minor_locator(LogLocator(subs='all', numticks=20)) -ax2.yaxis.set_minor_locator(LogLocator(subs='all', numticks=20)) -ax1.yaxis.set_minor_formatter(NullFormatter()) -ax2.yaxis.set_minor_formatter(NullFormatter()) -plt.grid(which='both', c='w') +lm_diff = np.abs((lm_aarhus_con - lm_empymod_con)) +ax2.loglog(lm_off_time, lm_diff, "m.", label=r"$|\Delta_\mathrm{LM}|$") +hm_diff = np.abs((hm_aarhus_con - hm_empymod_con)) +ax2.loglog(hm_off_time, hm_diff, "b.", label=r"$|\Delta_\mathrm{HM}|$") + +# Legend +ax2.legend(title="Difference", loc=3) + +# Labels and Settings +for ax in axs: + ax.set_ylabel(r"$\mathrm{d}\mathrm{B}_\mathrm{z}\,/\,\mathrm{d}t$") + ax.set_xlabel("Time(s)") + ax.grid(True, which="both", axis="both") + ax.yaxis.get_minor_locator().numticks = 30 + ############################################################################### empymod.Report() From c874ac45e395d66bbc28fc4cce04c080d71a2d77 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 00:04:54 +0100 Subject: [PATCH 02/14] Remove martera.eu, vanished --- CHANGELOG.rst | 12 +++++++++--- CREDITS.rst | 4 ++-- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f239feb..562566c 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -15,9 +15,15 @@ v2.5.x - Modelling routines: - - Arbitrary waveforms: Time-domain modelling can new be done not only for - impulse (``signal=0``), step-on (``signal=1``), and step-off - (``signal=-1``), but for any arbitrary piecewise linear waveform. + - ``bipole``: + + - Arbitrary waveforms: Time-domain modelling can new be done not only for + impulse (``signal=0``), step-on (``signal=1``), and step-off + (``signal=-1``), but for any arbitrary piecewise linear waveform. + + - Bandpass filters can now be applied to the frequency domain result + through user-provided functions. + - Merged ``loop`` into ``bipole``; there is no need for a special routine. Simply use ``bipole`` with ``msrc='b'`` and ``mrec=True`` for the same effect. As a result, ``empymod.utils.check_bipole`` changed its signature. diff --git a/CREDITS.rst b/CREDITS.rst index 5581e65..64916b6 100644 --- a/CREDITS.rst +++ b/CREDITS.rst @@ -16,8 +16,8 @@ namely: Aleksander Mousatov. - 2018-2021: `Delft University of Technology `_; through the *Gitaro.JIM* project (till 05/2021, emg3d v1.0.0), funded by - `MarTERA `_ as part of Horizon 2020, a funding scheme - of the European Research Area. + MarTERA as part of Horizon 2020, a funding scheme of the European Research + Area. - 2021-2022: `TERRASYS Geophysics GmbH & Co. KG `_. - 2021-2024: `Delft University of Technology `_ through From 9136cd6ad521f951324ee3062f03bf151e75e583 Mon Sep 17 00:00:00 2001 From: prisae Date: Thu, 29 Jan 2026 09:03:12 +0100 Subject: [PATCH 03/14] Update tudelft links to resolver --- docs/manual/references.rst | 2 +- examples/time_domain/step_and_impulse.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/manual/references.rst b/docs/manual/references.rst index 2ce06b6..9fe0455 100644 --- a/docs/manual/references.rst +++ b/docs/manual/references.rst @@ -20,7 +20,7 @@ References .. [Ghos70] Ghosh, D. P., 1970, The application of linear filter theory to the direct interpretation of geoelectrical resistivity measurements: Ph.D. Thesis, TU Delft; UUID: `88a568bb-ebee-4d7b-92df-6639b42da2b2 - `_. + `_. .. [GuSi97] Guptasarma, D., and B. Singh, 1997, New digital linear filters for Hankel J0 and J1 transforms: Geophysical Prospecting, 45, 745--762; DOI: `10.1046/j.1365-2478.1997.500292.x diff --git a/examples/time_domain/step_and_impulse.py b/examples/time_domain/step_and_impulse.py index 88853e2..56bb4cb 100644 --- a/examples/time_domain/step_and_impulse.py +++ b/examples/time_domain/step_and_impulse.py @@ -59,7 +59,7 @@ # Equations 3.2 and 3.3 in Werthmüller, D., 2009, Inversion of multi-transient # EM data from anisotropic media: M.S. thesis, TU Delft, ETH Zürich, RWTH # Aachen; UUID: `f4b071c1-8e55-4ec5-86c6-a2d54c3eda5a -# `_. +# `_. # # Analytical functions # ~~~~~~~~~~~~~~~~~~~~ From b1d2243926c45a742ef971ca5986ca7758ff19c9 Mon Sep 17 00:00:00 2001 From: prisae Date: Thu, 29 Jan 2026 09:07:21 +0100 Subject: [PATCH 04/14] rm marineemlab link; often fails, doi is enough --- docs/manual/references.rst | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/manual/references.rst b/docs/manual/references.rst index 9fe0455..98aabfd 100644 --- a/docs/manual/references.rst +++ b/docs/manual/references.rst @@ -42,9 +42,7 @@ References .. [Key09] Key, K., 2009, 1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers: Geophysics, 74(2), F9--F20; DOI: `10.1190/1.3058434 - `_. Software: - `marineemlab.ucsd.edu/Projects/Occam/1DCSEM - `_. + `_. .. [Key12] Key, K., 2012, Is the fast Hankel transform faster than quadrature?: Geophysics, 77(3), F21--F30; DOI: `10.1190/geo2011-0237.1 `_; Software: From 23d20f2375060829677f3c2c58dbbb4c167ea924 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 20:27:09 +0100 Subject: [PATCH 05/14] move waveform to tem --- empymod/model.py | 36 ++++++++++++++++++++++-------------- empymod/utils.py | 35 +++++++++++++++++++++++------------ 2 files changed, 45 insertions(+), 26 deletions(-) diff --git a/empymod/model.py b/empymod/model.py index 14c20d5..505a03a 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -427,8 +427,8 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, if signal is None: freq = freqtime else: - outtime = check_time(freqtime, signal, ft, ftarg, verb, True) - time, freq, ft, ftarg, signal, waveform = outtime + outtime = check_time(freqtime, signal, ft, ftarg, verb, new=True) + time, freq, ft, ftarg, signal = outtime # Check layer parameters model = check_model(depth, res, aniso, epermH, epermV, mpermH, mpermV, @@ -617,17 +617,6 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, # In case of QWE/QUAD, print Warning if not converged conv_warning(conv, ftarg, 'Fourier', verb) - # Apply waveform - if waveform: - map_time = waveform['map_time'][:, :, :] - wave_weight = waveform['wave_weight'][:, :, None] - wave_index = waveform['wave_index'][:, :, None] - gauss_weight = waveform['gauss_weight'][None, :, None, None] - - # Convert to time domain and apply waveform - EM = np.sum(np.sum(EM[map_time] * gauss_weight, axis=1) * - wave_weight * wave_index, axis=1) - # Reshape for number of sources EM = EM.reshape((-1, nrec, nsrc), order='F') @@ -1643,6 +1632,12 @@ def tem(fEM, off, freq, time, signal, ft, ftarg, conv=True): it can speed-up the calculation by omitting input-checks. """ + if isinstance(signal, dict): + waveform = signal + signal = waveform['signal'] + else: + waveform = False + # 1. Scale frequencies if switch-on/off response # Step function for causal times is like a unit fct, therefore an impulse # in frequency domain @@ -1660,4 +1655,17 @@ def tem(fEM, off, freq, time, signal, ft, ftarg, conv=True): tEM[:, i] += out[0] conv *= out[1] - return tEM*2/np.pi, conv # Scaling from Fourier transform + tEM *= 2/np.pi # Scaling from Fourier transform + + # Apply waveform + if waveform: + map_time = waveform['map_time'][:, :, :] + wave_weight = waveform['wave_weight'][:, :, None] + wave_index = waveform['wave_index'][:, :, None] + gauss_weight = waveform['gauss_weight'][None, :, None, None] + + # Convert to time domain and apply waveform + tEM = np.sum(np.sum(tEM[map_time] * gauss_weight, axis=1) * + wave_weight * wave_index, axis=1) + + return tEM, conv diff --git a/empymod/utils.py b/empymod/utils.py index cfcd6e1..99d12bc 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -1005,13 +1005,14 @@ def check_time(time, signal, ft, ftarg, verb, new=False): """ - if new and isinstance(signal, dict): - time, signal, waveform = check_waveform(time, **signal) - else: - waveform = None - # Check time and input signal - time = check_time_only(time, signal, verb) + out = check_time_only(time, signal, verb, new) + if new: + time, waveform = out + if isinstance(waveform, dict): + signal = waveform['signal'] + else: + time = out # Ensure ft is all lowercase ft = ft.lower() @@ -1256,12 +1257,12 @@ def check_time(time, signal, ft, ftarg, verb, new=False): # Make backwards compatible with old signature possibility? if new: - return time, freq, ft, targ, signal, waveform + return time, freq, ft, targ, waveform else: return time, freq, ft, targ -def check_time_only(time, signal, verb): +def check_time_only(time, signal, verb, new=False): r"""Check time and signal parameters. This check-function is called from one of the modelling routines in @@ -1294,9 +1295,15 @@ def check_time_only(time, signal, verb): global _min_time # Check input signal - if int(signal) not in [-1, 0, 1]: - raise ValueError(" must be one of: [None, -1, 0, 1]; " + if isinstance(signal, dict) and new: + time, signal = check_waveform(time, **signal) + # TODO print waveform info + elif int(signal) not in [-1, 0, 1]: + raise ValueError(" must be one of: [None, -1, 0, 1, dict]; " f" provided: {signal}") + else: + pass + # TODO print signal info # Check time time = _check_var(time, float, 1, 'time') @@ -1307,7 +1314,10 @@ def check_time_only(time, signal, verb): if verb > 2: _prnt_min_max_val(time, " time [s] : ", verb) - return time + if new: + return time, signal + else: + return time def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): @@ -1365,13 +1375,14 @@ def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): comp_time_flat, map_time = np.unique(comp_time, return_inverse=True) waveform = { + 'signal': signal, 'map_time': map_time, 'wave_weight': wave_weight, 'wave_index': wave_index, 'gauss_weight': gauss_weight, } - return comp_time_flat, signal, waveform + return comp_time_flat, waveform def check_solution(solution, signal, ab, msrc, mrec): From 161550dfcc7ab61ed727b45ba35a892f4015e8a5 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 21:58:57 +0100 Subject: [PATCH 06/14] Pass waveform fct --- empymod/model.py | 48 ++++++++++++++++++++++++++++++------------------ empymod/utils.py | 41 ++++++++++++++++++++--------------------- 2 files changed, 50 insertions(+), 39 deletions(-) diff --git a/empymod/model.py b/empymod/model.py index 505a03a..3e39333 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -137,13 +137,26 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, freqtime : array_like Frequencies f (Hz) if ``signal==None``, else times t (s); (f, t > 0). - signal : {None, 0, 1, -1}, default: None + signal : {None, -1, 0, 1, dict}, default: None Source signal: - None: Frequency-domain response - -1 : Switch-off time-domain response - 0 : Impulse time-domain response - +1 : Switch-on time-domain response + - dict : Arbitrary waveform + + For an arbitrary waveform, the dictionary must contain the following + keyword-value pairs: + + - nodes : array_like + Nodes of the waveform. + - amplitudes : array_like + Amplitudes (current) of the waveform. + - signal : {-1, 0, 1}, default: 0 + Signal that is convolved with the waveform. + - nquad : int, default: 3 + Number of quadrature points for the waveform segments. aniso : array_like, default: ones Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res. @@ -325,8 +338,16 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, If True, the output is squeezed. If False, the output will always be of ``ndim=3``, (nfreqtime, nrec, nsrc). - bandbass : dict or None, default: None - TODO + bandpass : {dict, None}, default: None + A dictionary containing any function that is applied to the + frequency-domain result. The signature of the function must be + ``func(inp, p_dict)``, where ``inp`` is the dictionary you provide, and + ``p_dict`` is a dictionary that contains all parameters so far computed + in empymod ``[locals()]``. Any change to the frequency domain result + must be done in-place, and the function does not return anything. Refer + to the time-domain loop examples in the gallery. The dictionary must + contain at least the keyword ``'func'``, containing the actual + function, but can contain any other parameters too. Returns @@ -407,7 +428,6 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, Out[2]: (1.6880934577857306e-10-3.083031298956568e-10j) """ - # Get kwargs with defaults. out = get_kwargs( [ @@ -619,8 +639,6 @@ def bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, # Reshape for number of sources EM = EM.reshape((-1, nrec, nsrc), order='F') - - # Squeeze empty dimensions if squeeze: EM = np.squeeze(EM) @@ -1632,11 +1650,12 @@ def tem(fEM, off, freq, time, signal, ft, ftarg, conv=True): it can speed-up the calculation by omitting input-checks. """ + # Check if standard signal or arbitrary waveform if isinstance(signal, dict): - waveform = signal - signal = waveform['signal'] + apply_waveform = signal['apply_waveform'] + signal = signal['signal'] else: - waveform = False + apply_waveform = False # 1. Scale frequencies if switch-on/off response # Step function for causal times is like a unit fct, therefore an impulse @@ -1658,14 +1677,7 @@ def tem(fEM, off, freq, time, signal, ft, ftarg, conv=True): tEM *= 2/np.pi # Scaling from Fourier transform # Apply waveform - if waveform: - map_time = waveform['map_time'][:, :, :] - wave_weight = waveform['wave_weight'][:, :, None] - wave_index = waveform['wave_index'][:, :, None] - gauss_weight = waveform['gauss_weight'][None, :, None, None] - - # Convert to time domain and apply waveform - tEM = np.sum(np.sum(tEM[map_time] * gauss_weight, axis=1) * - wave_weight * wave_index, axis=1) + if apply_waveform: + tEM = apply_waveform(tEM) return tEM, conv diff --git a/empymod/utils.py b/empymod/utils.py index 99d12bc..4d5d8f6 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -1336,8 +1336,7 @@ def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): # Pre-allocate comp_time = np.zeros((time.size, nquad, dIdt.size)) - wave_weight = np.zeros((time.size, dIdt.size)) - wave_index = np.zeros((time.size, dIdt.size), dtype=bool) + segment_weight = np.zeros((time.size, dIdt.size)) # Loop over wave segments. wi = 0 @@ -1348,41 +1347,41 @@ def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): continue # If wanted time is before a wave element, ignore it. - wave_ind = nodes[i] < time - if wave_ind.sum() == 0: + ind = nodes[i] < time + if ind.sum() == 0: continue # Start and end for this wave-segment for all times. - ta = time[wave_ind] - nodes[i] - tb = time[wave_ind] - nodes[i+1] + ta = time[ind] - nodes[i] + tb = time[ind] - nodes[i+1] # If wanted time is within a wave element, we cut the element. - tb[nodes[i+1] > time[wave_ind]] = 0.0 # Cut elements + tb[nodes[i+1] > time[ind]] = 0.0 # Cut elements # Gauss-Legendre for this wave segment. See # https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval # for the change of interval, which makes this a bit more complex. - comp_time[:, :, wi] = np.outer((tb-ta)/2, g_x)+(ta+tb)[:, None]/2 + comp_time[ind, :, wi] = np.outer((tb-ta)/2, g_x)+(ta+tb)[:, None]/2 - wave_weight[:, wi] = (tb-ta)/2*cdIdt - wave_index[:, wi] = wave_ind + segment_weight[ind, wi] = (tb-ta)/2*cdIdt wi += 1 comp_time = comp_time[:, :, :wi] - wave_weight = wave_weight[:, :wi] - wave_index = wave_index[:, :wi] + segment_weight = segment_weight[:, :wi] comp_time_flat, map_time = np.unique(comp_time, return_inverse=True) - waveform = { - 'signal': signal, - 'map_time': map_time, - 'wave_weight': wave_weight, - 'wave_index': wave_index, - 'gauss_weight': gauss_weight, - } - - return comp_time_flat, waveform + + def apply_waveform(resp): + """Waveform function.""" + + # Gauss quadrature of each segment + resp = np.sum(resp[map_time]*gauss_weight[None, :, None, None], axis=1) + + # Sum over waveform elements + return np.sum(resp*segment_weight[:, :, None], axis=1) + + return comp_time_flat, {'signal': signal, 'apply_waveform': apply_waveform} def check_solution(solution, signal, ab, msrc, mrec): From 93325d0376c4e85768e61ef64ba32c35914e171b Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 23:14:57 +0100 Subject: [PATCH 07/14] Clean up & document utils.py --- CHANGELOG.rst | 34 ++++++++----- empymod/model.py | 7 +-- empymod/utils.py | 116 ++++++++++++++++++++++++++++++++++++-------- tests/test_utils.py | 84 ++++++++++++++++++++------------ 4 files changed, 175 insertions(+), 66 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 562566c..c24040c 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -17,21 +17,31 @@ v2.5.x - ``bipole``: - - Arbitrary waveforms: Time-domain modelling can new be done not only for - impulse (``signal=0``), step-on (``signal=1``), and step-off - (``signal=-1``), but for any arbitrary piecewise linear waveform. - - - Bandpass filters can now be applied to the frequency domain result - through user-provided functions. + - Arbitrary waveforms: Time-domain modelling can new be done for any + arbitrary piecewise linear waveform, in addition to impulse + (``signal=0``), step-on (``signal=1``), and step-off (``signal=-1``) + responses. + - User-defined bandpass filters can now be applied to the frequency domain + result through user-provided functions. - Merged ``loop`` into ``bipole``; there is no need for a special routine. Simply use ``bipole`` with ``msrc='b'`` and ``mrec=True`` for the same - effect. As a result, ``empymod.utils.check_bipole`` changed its signature. - For ``bipole``, it new prints the source and receiver type (electric field, - magnetic field, magnetic flux, or electric current). Both, - ``empymod.model.loop`` and ``empymod.utils.check_bipole``'s old signature - are deprecated and will be removed in v3. + effect. + +- New prints (if verbose): + + - The source/receiver types are new printed. + - The signal/waveform is new printed. + +- New function ``empymod.utils.check_waveform``, to check the waveform. + +- New deprecations (will be removed in v3): + - ``empymod.model.loop`` will be removed in v3. + - ``empymod.utils.check_bipole``: Now returns also the src/rec field and + type. + - ``empymod.utils.check_time``: Now also returns the signal. + - ``empymod.utils.check_time_only``: Now also returns the signal. v2.5.4: Bugfix ext. fcts with z+ up @@ -152,7 +162,7 @@ The code is now compatible with NumPy v2. - Gallery Update Part I: - - Update for Jupyterlab (ipympl/widget) + - Update for Jupyterlab (ipympl/widget) - Replaced implicit by explicit pyplots - Use by default a positiv z-upwards coordinate system - Part I: frequency domain; reproducing; published diff --git a/empymod/model.py b/empymod/model.py index 3e39333..4d1751a 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -820,7 +820,8 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, # Check times and Fourier Transform arguments, get required frequencies # (freq = freqtime if `signal=None`) if signal is not None: - time, freq, ft, ftarg = check_time(freqtime, signal, ft, ftarg, verb) + out = check_time(freqtime, signal, ft, ftarg, verb, True) + time, freq, ft, ftarg, _ = out else: freq = freqtime @@ -1084,7 +1085,7 @@ def analytical(src, rec, res, freqtime, solution='fs', signal=None, ab=11, # Check times or frequencies if signal is not None: - freqtime = check_time_only(freqtime, signal, verb) + freqtime, _ = check_time_only(freqtime, signal, verb, True) # Check layer parameters model = check_model([], res, aniso, epermH, epermV, mpermH, mpermV, True, @@ -1218,7 +1219,7 @@ def gpr(src, rec, depth, res, freqtime, cf, gain=None, ab=11, aniso=None, # === 1. CHECK TIME ============ # Check times and Fourier Transform arguments, get required frequencies - time, freq, ft, ftarg = check_time(freqtime, 0, ft, ftarg, verb) + time, freq, ft, ftarg, _ = check_time(freqtime, 0, ft, ftarg, verb, True) # === 2. CALL DIPOLE ============ diff --git a/empymod/utils.py b/empymod/utils.py index 4d5d8f6..48649c0 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -973,13 +973,13 @@ def check_time(time, signal, ft, ftarg, verb, new=False): time : array_like Times t (s). - signal : {None, 0, 1, -1} + signal : {-1, 0, 1, dict} Source signal: - - None: Frequency-domain response - -1 : Switch-off time-domain response - 0 : Impulse time-domain response - +1 : Switch-on time-domain response + - dict : Arbitrary waveform ft : {'dlf', 'qwe', 'fftlog', 'fft'} Flag for Fourier transform. @@ -990,6 +990,9 @@ def check_time(time, signal, ft, ftarg, verb, new=False): verb : {0, 1, 2, 3, 4} Level of verbosity. + new : bool, default=False + For backwards compatibility. + Returns ------- @@ -1003,6 +1006,10 @@ def check_time(time, signal, ft, ftarg, verb, new=False): Checked if valid and set to defaults if not provided, checked with signal. + signal + Signal, checked for input and with waveform function, if applicable. + Only returned if new=True. + """ # Check time and input signal @@ -1255,10 +1262,16 @@ def check_time(time, signal, ft, ftarg, verb, new=False): if args and verb > 0: print(f"* WARNING :: Unknown ftarg {args} for method '{ft}'") - # Make backwards compatible with old signature possibility? + # Backwards compatible return. if new: return time, freq, ft, targ, waveform else: + msg = ( + "The signature of `empymod.utils.check_time` changed from " + "returning `time, freq, ft, targ` to returning `time, freq, ft, " + "targ, waveform`. The old signature will be removed in v3.0." + ) + warnings.warn(msg, DeprecationWarning) return time, freq, ft, targ @@ -1274,37 +1287,33 @@ def check_time_only(time, signal, verb, new=False): time : array_like Times t (s). - signal : {None, 0, 1, -1} + signal : {-1, 0, 1, dict} Source signal: - - None: Frequency-domain response - -1 : Switch-off time-domain response - 0 : Impulse time-domain response - +1 : Switch-on time-domain response + - dict : Arbitrary waveform verb : {0, 1, 2, 3, 4} Level of verbosity. + new : bool, default=False + For backwards compatibility. + Returns ------- time : float Time, checked for size and assured min_time. + signal + Signal, checked for input and with waveform function, if applicable. + Only returned if new=True. + """ global _min_time - # Check input signal - if isinstance(signal, dict) and new: - time, signal = check_waveform(time, **signal) - # TODO print waveform info - elif int(signal) not in [-1, 0, 1]: - raise ValueError(" must be one of: [None, -1, 0, 1, dict]; " - f" provided: {signal}") - else: - pass - # TODO print signal info - # Check time time = _check_var(time, float, 1, 'time') @@ -1314,16 +1323,75 @@ def check_time_only(time, signal, verb, new=False): if verb > 2: _prnt_min_max_val(time, " time [s] : ", verb) + # Check input signal + if isinstance(signal, dict) and new: + time, signal = check_waveform(time, **signal, verb=verb) + if verb > 2: + print(f" signal : {signal['signal']}") + + elif int(signal) not in [-1, 0, 1]: + raise ValueError(" must be one of: [None, -1, 0, 1, dict]; " + f" provided: {signal}") + elif verb > 2: + print(f" signal : {signal}") + + # Backwards compatible return. if new: return time, signal else: + msg = ( + "The signature of `empymod.utils.check_time_only` changed from " + "returning `time` to returning `time, signal`. " + "The old signature will be removed in v3.0." + ) + warnings.warn(msg, DeprecationWarning) return time -def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): +def check_waveform(time, nodes, amplitudes, verb, signal=0, nquad=3): + r"""Check waveform dict and create waveform function. - # The waveform function (here and in model.bipole) is based on work - # from Kerry Key. + This check-function is called from :mod:`empymod.utils.check_time_only`. + + The waveform function is based on work from Kerry Key (DIPOLE1D). + + + Parameters + ---------- + time : array_like + Times t (s). + + nodes : array_like + Waveform nodes, for which amplitudes are provided. + + amplitudes : array_like + Waveform amplitudes corresponding to the nodes. + + verb : {0, 1, 2, 3, 4} + Level of verbosity. + + signal : {-1, 0, 1}, default: 0 + Source signal to use for convolution: + + - -1 : Switch-off time-domain response + - 0 : Impulse time-domain response + - +1 : Switch-on time-domain response + + nquad : int, default : 3 + Quadrature points for wave segments. + + + Returns + ------- + time : float + Time required to compute for the waveform. + + signal + Signal dictionary, containing the signal itself and the waveform + function. + + """ + global _min_time # Waveform segments. dt = np.diff(nodes) @@ -1356,7 +1424,7 @@ def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): tb = time[ind] - nodes[i+1] # If wanted time is within a wave element, we cut the element. - tb[nodes[i+1] > time[ind]] = 0.0 # Cut elements + tb[nodes[i+1] > time[ind]] = 0.0 # Gauss-Legendre for this wave segment. See # https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval @@ -1367,10 +1435,13 @@ def check_waveform(time, nodes, amplitudes, signal=0, nquad=3): wi += 1 + # Cut arrays to number of relevant segments. comp_time = comp_time[:, :, :wi] segment_weight = segment_weight[:, :wi] + # Get unique times to compute, with reverse indices. comp_time_flat, map_time = np.unique(comp_time, return_inverse=True) + comp_time_flat[comp_time_flat < _min_time] = _min_time def apply_waveform(resp): """Waveform function.""" @@ -1381,6 +1452,11 @@ def apply_waveform(resp): # Sum over waveform elements return np.sum(resp*segment_weight[:, :, None], axis=1) + if verb > 2: + _prnt_min_max_val(comp_time_flat, " time comp. [s] : ", verb) + print(f" wave nodes [s] : {_strvar(nodes)}") + print(f" wave ampl. [-] : {_strvar(amplitudes)}") + return comp_time_flat, {'signal': signal, 'apply_waveform': apply_waveform} diff --git a/tests/test_utils.py b/tests/test_utils.py index dc9ad4b..403fe7e 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -559,7 +559,8 @@ def test_check_time(capsys): # # DLF # # # verbose - _, f, ft, ftarg = utils.check_time(time, 0, 'dlf', {}, 4) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, ft, ftarg = utils.check_time(time, 0, 'dlf', {}, 4) out, _ = capsys.readouterr() assert " time [s] : 3" in out assert " Fourier : DLF (Sine-Filter)" in out @@ -579,12 +580,14 @@ def test_check_time(capsys): assert ftarg['kind'] == 'sin' # filter-string and unknown parameter - _, f, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, _, ftarg = utils.check_time( time, -1, 'dlf', {'dlf': 'key_201_2012', 'kind': 'cos', 'notused': 1}, 4) out, _ = capsys.readouterr() outstr = " time [s] : 3\n" + outstr = " signal : -1\n" outstr += " Fourier : DLF (Cosine-Filter)\n > Filter" assert outstr in out assert "WARNING :: Unknown ftarg {'notused': 1} for method 'dlf'" in out @@ -597,7 +600,8 @@ def test_check_time(capsys): assert ftarg['kind'] == 'cos' # filter instance - _, _, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time( time, 1, 'dlf', {'dlf': filters.Fourier().key_201_2012, 'kind': 'sin'}, 0) assert ftarg['dlf'].name == filters.Fourier().key_201_2012.name @@ -606,7 +610,8 @@ def test_check_time(capsys): # pts_per_dec out, _ = capsys.readouterr() # clear buffer - _, _, _, ftarg = utils.check_time(time, 0, 'dlf', {'pts_per_dec': 30}, 4) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time(time, 0, 'dlf', {'pts_per_dec': 30}, 4) assert ftarg['dlf'].name == filters.Fourier().key_201_2012.name assert ftarg['pts_per_dec'] == 30 assert ftarg['kind'] == 'sin' @@ -614,7 +619,8 @@ def test_check_time(capsys): assert " > DLF type : Splined, 30.0 pts/dec" in out # filter-string and pts_per_dec - _, _, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time( time, 0, 'dlf', {'dlf': 'key_81_2009', 'pts_per_dec': -1, 'kind': 'cos'}, 4) out, _ = capsys.readouterr() @@ -624,7 +630,8 @@ def test_check_time(capsys): assert ftarg['kind'] == 'cos' # pts_per_dec - _, freq, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, freq, _, ftarg = utils.check_time( time, 0, 'dlf', {'pts_per_dec': 0, 'kind': 'sin'}, 4) out, _ = capsys.readouterr() assert " > DLF type : Standard" in out @@ -633,7 +640,8 @@ def test_check_time(capsys): assert_allclose(np.ravel(f_base/(2*np.pi*time[:, None])), freq) # filter-string and pts_per_dec - _, _, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time( time, 0, 'dlf', {'dlf': 'key_81_2009', 'pts_per_dec': 50, 'kind': 'cos'}, 0) assert ftarg['dlf'].name == filters.Fourier().key_81_2009.name @@ -641,7 +649,8 @@ def test_check_time(capsys): assert ftarg['kind'] == 'cos' # just kind - _, f, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, _, ftarg = utils.check_time( time, 0, 'dlf', {'kind': 'sin'}, 0) assert ftarg['pts_per_dec'] == -1 assert_allclose(f[:9], f1) @@ -650,17 +659,20 @@ def test_check_time(capsys): # Assert it can be called repetitively _, _ = capsys.readouterr() - _, _, _, ftarg = utils.check_time(time, 0, 'dlf', {}, 1) - _, _, _, _ = utils.check_time(time, 0, 'dlf', ftarg, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time(time, 0, 'dlf', {}, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, _ = utils.check_time(time, 0, 'dlf', ftarg, 1) out, _ = capsys.readouterr() assert out == "" # # QWE # # # verbose - _, f, ft, ftarg = utils.check_time(time, 0, 'qwe', {}, 4) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, ft, ftarg = utils.check_time(time, 0, 'qwe', {}, 4) out, _ = capsys.readouterr() outstr = " Fourier : Quadrature-with-Extrapolation\n > rtol" - assert out[24:87] == outstr + assert out[48:111] == outstr assert ft == 'qwe' assert ftarg['rtol'] == 1e-8 assert ftarg['atol'] == 1e-20 @@ -683,7 +695,8 @@ def test_check_time(capsys): assert ftarg['sincos'] is np.sin # only limit - _, _, _, ftarg = utils.check_time(time, 1, 'qwe', {'limit': 30}, 0) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time(time, 1, 'qwe', {'limit': 30}, 0) assert ftarg['rtol'] == 1e-8 assert ftarg['atol'] == 1e-20 assert ftarg['nquad'] == 21 @@ -696,7 +709,8 @@ def test_check_time(capsys): assert ftarg['sincos'] is np.sin # all arguments - _, _, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time( time, -1, 'qwe', {'rtol': 1e-3, 'atol': 1e-4, 'nquad': 31, 'maxint': 20, 'pts_per_dec': 30, 'diff_quad': 200, 'a': 0.01, 'b': 0.2, @@ -719,14 +733,17 @@ def test_check_time(capsys): # Assert it can be called repetitively _, _ = capsys.readouterr() - _, _, _, ftarg = utils.check_time(time, 0, 'qwe', {}, 1) - _, _, _, _ = utils.check_time(time, 0, 'qwe', ftarg, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time(time, 0, 'qwe', {}, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, _ = utils.check_time(time, 0, 'qwe', ftarg, 1) out, _ = capsys.readouterr() assert out == "" # # FFTLog # # # verbose - _, f, ft, ftarg = utils.check_time(time, 0, 'fftlog', {}, 4) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, ft, ftarg = utils.check_time(time, 0, 'fftlog', {}, 4) out, _ = capsys.readouterr() outstr = " Fourier : FFTLog\n > pts_per_dec" assert outstr in out @@ -759,7 +776,8 @@ def test_check_time(capsys): assert_allclose(f, fres, rtol=1e-5) # Several parameters - _, _, _, ftarg = utils.check_time( + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time( time, -1, 'fftlog', {'pts_per_dec': 10, 'add_dec': [-3, 4], 'q': 2}, 0) assert ftarg['pts_per_dec'] == 10 @@ -772,14 +790,17 @@ def test_check_time(capsys): # Assert it can be called repetitively _, _ = capsys.readouterr() - _, _, _, ftarg = utils.check_time(time, 0, 'fftlog', {}, 1) - _, _, _, _ = utils.check_time(time, 0, 'fftlog', ftarg, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, ftarg = utils.check_time(time, 0, 'fftlog', {}, 1) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, _, _, _ = utils.check_time(time, 0, 'fftlog', ftarg, 1) out, _ = capsys.readouterr() assert out == "" # # FFT # # # verbose - _, f, ft, ftarg = utils.check_time(time, 0, 'fft', {}, 4) + with pytest.warns(DeprecationWarning, match='in v3.0.'): + _, f, ft, ftarg = utils.check_time(time, 0, 'fft', {}, 4) out, _ = capsys.readouterr() assert "Fourier : Fast Fourier Transform FFT\n > dfreq" in out assert " > pts_per_dec : (linear)" in out @@ -794,14 +815,15 @@ def test_check_time(capsys): assert_allclose(f[-5:], fres[-5:]) # Several parameters - _, _, _, ftarg = utils.check_time( - time, 0, 'fft', {'dfreq': 1e-3, 'nfreq': 2**15+1, 'ntot': 3}, 0) + _, _, _, ftarg, _ = utils.check_time( + time, 0, 'fft', {'dfreq': 1e-3, 'nfreq': 2**15+1, 'ntot': 3}, 0, True) assert ftarg['dfreq'] == 0.001 assert ftarg['nfreq'] == 2**15+1 assert ftarg['ntot'] == 2**16 # Several parameters; pts_per_dec - _, f, _, ftarg = utils.check_time(time, 0, 'fft', {'pts_per_dec': 5}, 3) + _, f, _, ftarg, _ = utils.check_time( + time, 0, 'fft', {'pts_per_dec': 5}, 3, True) out, _ = capsys.readouterr() assert " > pts_per_dec : 5" in out assert ftarg['dfreq'] == 0.002 @@ -819,8 +841,8 @@ def test_check_time(capsys): # Assert it can be called repetitively _, _ = capsys.readouterr() - _, _, _, ftarg = utils.check_time(time, 0, 'fft', {}, 1) - _, _, _, _ = utils.check_time(time, 0, 'fft', ftarg, 1) + _, _, _, ftarg, _ = utils.check_time(time, 0, 'fft', {}, 1, True) + _, _, _, _, _ = utils.check_time(time, 0, 'fft', ftarg, 1, True) out, _ = capsys.readouterr() assert out == "" @@ -828,26 +850,26 @@ def test_check_time(capsys): # minimum time _ = utils.check_time( - 0, 0, 'dlf', {'dlf': 'key_201_2012', 'kind': 'cos'}, 1) + 0, 0, 'dlf', {'dlf': 'key_201_2012', 'kind': 'cos'}, 1, True) out, _ = capsys.readouterr() assert out[:21] == "* WARNING :: Times < " # Signal != -1, 0, 1 with pytest.raises(ValueError, match=' must be one of:'): - utils.check_time(time, -2, 'dlf', {}, 0) + utils.check_time(time, -2, 'dlf', {}, 0, True) # ft != cos, sin, dlf, qwe, fftlog, with pytest.raises(ValueError, match=' must be one of:'): - utils.check_time(time, 0, 'bla', {}, 0) + utils.check_time(time, 0, 'bla', {}, 0, True) # filter missing attributes dlf = filters.DigitalFilter('test') with pytest.raises(AttributeError, match='DLF-filter is missing some att'): - utils.check_time(time, 0, 'dlf', {'dlf': dlf}, 1) + utils.check_time(time, 0, 'dlf', {'dlf': dlf}, 1, True) # filter with wrong kind with pytest.raises(ValueError, match="'kind' must be either 'sin' or"): - utils.check_time(time, 0, 'dlf', {'kind': 'wrongkind'}, 1) + utils.check_time(time, 0, 'dlf', {'kind': 'wrongkind'}, 1, True) def test_check_solution(capsys): From 4f37a9c90e17ff049d58c2be5fa2223cf74b2f42 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 23:24:39 +0100 Subject: [PATCH 08/14] Add arb. waveform & bandpass to dipole --- CHANGELOG.rst | 2 +- empymod/model.py | 41 ++++++++++++++++++++++++++++++++++++----- tests/test_utils.py | 3 ++- 3 files changed, 39 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index c24040c..01a71e2 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -15,7 +15,7 @@ v2.5.x - Modelling routines: - - ``bipole``: + - ``bipole`` and ``dipole``: - Arbitrary waveforms: Time-domain modelling can new be done for any arbitrary piecewise linear waveform, in addition to impulse diff --git a/empymod/model.py b/empymod/model.py index 4d1751a..81c0b17 100644 --- a/empymod/model.py +++ b/empymod/model.py @@ -714,13 +714,26 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, freqtime : array_like Frequencies f (Hz) if ``signal==None``, else times t (s); (f, t > 0). - signal : {None, 0, 1, -1}, default: None + signal : {None, -1, 0, 1, dict}, default: None Source signal: - None: Frequency-domain response - -1 : Switch-off time-domain response - 0 : Impulse time-domain response - +1 : Switch-on time-domain response + - dict : Arbitrary waveform + + For an arbitrary waveform, the dictionary must contain the following + keyword-value pairs: + + - nodes : array_like + Nodes of the waveform. + - amplitudes : array_like + Amplitudes (current) of the waveform. + - signal : {-1, 0, 1}, default: 0 + Signal that is convolved with the waveform. + - nquad : int, default: 3 + Number of quadrature points for the waveform segments. ab : int, default: 11 Source-receiver configuration. @@ -772,6 +785,17 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, If True, the output is squeezed. If False, the output will always be of ``ndim=3``, (nfreqtime, nrec, nsrc). + bandpass : {dict, None}, default: None + A dictionary containing any function that is applied to the + frequency-domain result. The signature of the function must be + ``func(inp, p_dict)``, where ``inp`` is the dictionary you provide, and + ``p_dict`` is a dictionary that contains all parameters so far computed + in empymod ``[locals()]``. Any change to the frequency domain result + must be done in-place, and the function does not return anything. Refer + to the time-domain loop examples in the gallery. The dictionary must + contain at least the keyword ``'func'``, containing the actual + function, but can contain any other parameters too. + Returns ------- @@ -807,10 +831,13 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, """ # Get kwargs with defaults. out = get_kwargs( - ['verb', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'squeeze'], - [2, 'dlf', {}, 'dlf', {}, False, None, True], kwargs, + [ + 'verb', 'ht', 'htarg', 'ft', 'ftarg', 'xdirect', 'loop', 'squeeze', + 'bandpass', + ], + [2, 'dlf', {}, 'dlf', {}, False, None, True, None], kwargs, ) - verb, ht, htarg, ft, ftarg, xdirect, loop, squeeze = out + verb, ht, htarg, ft, ftarg, xdirect, loop, squeeze, bandpass = out # === 1. LET'S START ============ t0 = printstartfinish(verb) @@ -821,7 +848,7 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, # (freq = freqtime if `signal=None`) if signal is not None: out = check_time(freqtime, signal, ft, ftarg, verb, True) - time, freq, ft, ftarg, _ = out + time, freq, ft, ftarg, signal = out else: freq = freqtime @@ -873,6 +900,10 @@ def dipole(src, rec, depth, res, freqtime, signal=None, ab=11, aniso=None, # In case of QWE/QUAD, print Warning if not converged conv_warning(conv, htarg, 'Hankel', verb) + # Apply user-provided bandpass filter + if isinstance(bandpass, dict) and 'func' in bandpass: + bandpass['func'](bandpass, locals()) + # Do f->t transform if required if signal is not None: EM, conv = tem(EM, off, freq, time, signal, ft, ftarg) diff --git a/tests/test_utils.py b/tests/test_utils.py index 403fe7e..268a3cf 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -611,7 +611,8 @@ def test_check_time(capsys): # pts_per_dec out, _ = capsys.readouterr() # clear buffer with pytest.warns(DeprecationWarning, match='in v3.0.'): - _, _, _, ftarg = utils.check_time(time, 0, 'dlf', {'pts_per_dec': 30}, 4) + _, _, _, ftarg = utils.check_time( + time, 0, 'dlf', {'pts_per_dec': 30}, 4) assert ftarg['dlf'].name == filters.Fourier().key_201_2012.name assert ftarg['pts_per_dec'] == 30 assert ftarg['kind'] == 'sin' From c60d0bb774f1ace843f57fd6d64023bd196849ff Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Thu, 29 Jan 2026 23:29:56 +0100 Subject: [PATCH 09/14] Update kwargs info --- empymod/utils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/empymod/utils.py b/empymod/utils.py index 48649c0..ce025f3 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -1979,15 +1979,15 @@ def get_kwargs(names, defaults, kwargs): - ALL functions: src, rec, res, aniso, epermH, epermV, mpermH, mpermV, verb - ONLY gpr: cf, gain - - ONLY bipole: msrc, srcpts, bandpass + - ONLY bipole: msrc, srcpts, mrec, recpts, strength + - ONLY bipole, dipole: bandpass - ONLY dipole_k: freq, wavenumber - ONLY analytical: solution - - ONLY bipole, loop: mrec, recpts, strength - - ONLY bipole, dipole, loop, gpr: ht, htarg, ft, ftarg, xdirect, loop - - ONLY bipole, dipole, loop, analytical: signal, squeeze + - ONLY bipole, dipole, gpr: ht, htarg, ft, ftarg, xdirect, loop + - ONLY bipole, dipole, analytical: signal, squeeze - ONLY dipole, analytical, gpr, dipole_k: ab - - ONLY bipole, dipole, loop, gpr, dipole_k: depth - - ONLY bipole, dipole, loop, analytical, gpr: freqtime + - ONLY bipole, dipole, gpr, dipole_k: depth + - ONLY bipole, dipole, analytical, gpr: freqtime Parameters ---------- From d1f910872fed66f668efd36bb0796e78640222bc Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Fri, 30 Jan 2026 15:20:33 +0100 Subject: [PATCH 10/14] Tests model --- tests/test_model.py | 88 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/tests/test_model.py b/tests/test_model.py index 41e9e18..2ed669c 100644 --- a/tests/test_model.py +++ b/tests/test_model.py @@ -693,6 +693,94 @@ def func_zeta(inp, pdict): assert_allclose(standard, outzeta) +def test_bandpass(): + def bandpass(inp, p_inp): + f = p_inp['freq'] + p_inp['EM'][f < 1] *= 0 + p_inp['EM'][f > 1] *= inp['fact'] + + mod = { + 'depth': [0], + 'res': [1.5, 5], + 'freqtime': np.array([0.1, 10]), + } + + # bipole() + bip = { + 'src': [0, 0, 0, 90, 0], + 'rec': [8000, 200, 300, 0, 0], + } + efield = bipole(**bip, **mod) + fact = 2 + bp = bipole(bandpass={'func': bandpass, 'fact': fact}, **bip, **mod) + assert_allclose(0+0j, bp[0]) + assert_allclose(fact*efield[1], bp[1]) + + # dipole() + dip = { + 'src': [0, 0, 0], + 'rec': [8000, 200, 300], + 'ab': 14, + } + efield = dipole(**dip, **mod) + fact = 2 + bp = dipole(bandpass={'func': bandpass, 'fact': fact}, **dip, **mod) + assert_allclose(0+0j, bp[0]) + assert_allclose(fact*efield[1], bp[1]) + + +def test_waveform(): + # Compare waveforms against signal=-1, 0, 1. + + time = np.logspace(-1, 1, 11) + + inp = { + 'src': [0, 0, 0.001, 0, 0], + 'rec': [6000, 0, 0.001, 0, 0], + 'depth': [0], + 'res': [2e14, 10], + 'aniso': [1, 2], + 'epermH': [0, 1], + 'xdirect': True, + 'freqtime': time, + 'ftarg': {'dlf': 'key_81_2009'}, + 'verb': 1, + } + + # Step-on : Factor -1 + signal = 1 + ramp = {'nodes': [-0.1e-6, 0.1e-6], 'amplitudes': [0, 1], 'signal': 1} + out1 = bipole(**inp, signal=signal) + out2 = bipole(**inp, signal=ramp) + assert_allclose(out1, -out2) + + # Step off + signal = -1 + ramp = {'nodes': [-0.1e-6, 0.1e-6], 'amplitudes': [1, 0], 'signal': -1} + out1 = bipole(**inp, signal=signal) + out2 = bipole(**inp, signal=ramp) + assert_allclose(out1, out2) + + # Impulse + signal = 0 + ramp1 = {'nodes': [-1e-5, -0.5e-6, 0, 0.5e-6, 1e-5], + 'amplitudes': [0, 0, 2e6, 0, 0], 'signal': 1} + ramp2 = {'nodes': [-1e-5, -0.5e-6, 0, 0.5e-6, 1e-5], + 'amplitudes': [0, 0, 2e6, 0, 0], 'signal': -1} + ramp3 = {'nodes': [-1e-6, 0], 'amplitudes': [1, 0], 'signal': 0} + ramp4 = {'nodes': [-1e-6, 0], 'amplitudes': [0, 1], 'signal': 0} + + out0 = bipole(**inp, signal=signal) + out1 = bipole(**inp, signal=ramp1) + out2 = bipole(**inp, signal=ramp2) + out3 = bipole(**inp, signal=ramp3) + out4 = bipole(**inp, signal=ramp4) + assert_allclose(out0, -out1, atol=2e-13) + assert_allclose(out0, out2, atol=2e-13) + assert_allclose(out0, out3, atol=2e-13) + assert_allclose(out0, -out4, atol=2e-13) + + def test_all_depths(): # Test RHS/LHS low-to-high/high-to-low src = [0, 0, 10] From dd03d5ec052b2de42e9cbe7e20a5602ba8b0cd5a Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Fri, 30 Jan 2026 16:22:15 +0100 Subject: [PATCH 11/14] Tests util --- empymod/utils.py | 5 ++++ tests/test_utils.py | 58 +++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 61 insertions(+), 2 deletions(-) diff --git a/empymod/utils.py b/empymod/utils.py index ce025f3..74cc87e 100644 --- a/empymod/utils.py +++ b/empymod/utils.py @@ -1443,6 +1443,11 @@ def check_waveform(time, nodes, amplitudes, verb, signal=0, nquad=3): comp_time_flat, map_time = np.unique(comp_time, return_inverse=True) comp_time_flat[comp_time_flat < _min_time] = _min_time + if comp_time_flat.size == 0: + raise ValueError( + "All wanted times are before provided waveform; aborting." + ) + def apply_waveform(resp): """Waveform function.""" diff --git a/tests/test_utils.py b/tests/test_utils.py index 268a3cf..8f205fc 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -816,14 +816,15 @@ def test_check_time(capsys): assert_allclose(f[-5:], fres[-5:]) # Several parameters - _, _, _, ftarg, _ = utils.check_time( + _, _, _, ftarg, sig = utils.check_time( time, 0, 'fft', {'dfreq': 1e-3, 'nfreq': 2**15+1, 'ntot': 3}, 0, True) assert ftarg['dfreq'] == 0.001 assert ftarg['nfreq'] == 2**15+1 assert ftarg['ntot'] == 2**16 + assert sig == 0 # Several parameters; pts_per_dec - _, f, _, ftarg, _ = utils.check_time( + _, f, _, ftarg, sig = utils.check_time( time, 0, 'fft', {'pts_per_dec': 5}, 3, True) out, _ = capsys.readouterr() assert " > pts_per_dec : 5" in out @@ -831,6 +832,7 @@ def test_check_time(capsys): assert ftarg['nfreq'] == 2048 assert ftarg['ntot'] == 2048 assert ftarg['pts_per_dec'] == 5 + assert sig == 0 outf = np.array([2.00000000e-03, 3.22098066e-03, 5.18735822e-03, 8.35419026e-03, 1.34543426e-02, 2.16680888e-02, 3.48962474e-02, 5.62000691e-02, 9.05096680e-02, @@ -873,6 +875,58 @@ def test_check_time(capsys): utils.check_time(time, 0, 'dlf', {'kind': 'wrongkind'}, 1, True) +def test_check_time_only(capsys): + time = np.array([1, 2, 3]) + + out, _ = capsys.readouterr() + waveform = { + 'nodes': [-1, -.8, -.2, 0], + 'amplitudes': [0, 1, 1, 0], + } + _, signal = utils.check_time_only(time, waveform, 3, True) + out, _ = capsys.readouterr() + assert 'signal : 0' in out + assert 'wave nodes [s] : -1 -0.8 -0.2 0' in out + assert 'wave ampl. [-] : 0 1 1 0' in out + assert signal['signal'] == 0 + + +def test_check_waveform(capsys): + time = np.array([.1, .2, .3]) + + out, _ = capsys.readouterr() + waveform = { + 'nodes': [-1, -.8, -.2, .15], + 'amplitudes': [0, 1, 1, 0], + } + + otime, signal = utils.check_waveform(time, **waveform, verb=3) + out, _ = capsys.readouterr() + assert 'time comp. [s] : 0.0338105 - 1.27746 : 18' in out + assert signal['signal'] == 0 + tEM = np.ones((otime.size, 1)) + EM = signal['apply_waveform'](tEM) + assert EM.shape == (time.size, 1) + print(EM[0][0]) + assert_allclose(EM, [[-0.14285714285714313], [0], [0]]) # Only status quo + + # empty comp times + with pytest.raises(ValueError, match='provided waveform; aborting'): + time = np.array([1]) + waveform = {'nodes': [2, 3], 'amplitudes': [0, 1], } + _ = utils.check_waveform(time, **waveform, verb=3) + + # some times before + time = np.array([.05, .1, .2]) + waveform = { + 'nodes': [.15, 1], + 'amplitudes': [0, 1], + } + otime, signal = utils.check_waveform(time, **waveform, verb=3) + # Only status quo + assert_allclose(otime, [1e-20, 5.635083e-3, 2.5e-2, 4.436492e-2]) + + def test_check_solution(capsys): # wrong solution with pytest.raises(ValueError, match='Solution must be one of'): From c6293bbd54ec1a688b2077f9f24b8f56bf601cea Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Fri, 30 Jan 2026 16:58:21 +0100 Subject: [PATCH 12/14] Bump minimum requirements --- .github/workflows/linux.yml | 14 +++++++------- CHANGELOG.rst | 6 ++++++ docs/manual/installation.rst | 4 ++-- pyproject.toml | 6 +++--- 4 files changed, 18 insertions(+), 12 deletions(-) diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index f6970a6..3a4670b 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -39,37 +39,37 @@ jobs: - python-version: "3.11" name: minimal os: ubuntu - conda: "numba 'scipy=1.10' scooby libdlf" + conda: "numba 'numpy==2.0.0' scipy scooby libdlf" test: "" - python-version: "3.12" name: plain os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf" test: "" - python-version: "3.12" name: full os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf matplotlib pytest-mpl" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf matplotlib pytest-mpl" test: "--mpl" - python-version: "3.13" name: plain os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf" test: "" - python-version: "3.13" name: full os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf matplotlib pytest-mpl" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf matplotlib pytest-mpl" test: "--mpl" - python-version: "3.14" name: plain os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf" test: "" - python-version: "3.14" name: full os: ubuntu - conda: "numba 'scipy>=1.10' scooby libdlf matplotlib pytest-mpl" + conda: "numba 'numpy>=2.0.0' scipy scooby libdlf matplotlib pytest-mpl" test: "--mpl" env: diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 01a71e2..b642c43 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -43,6 +43,12 @@ v2.5.x - ``empymod.utils.check_time``: Now also returns the signal. - ``empymod.utils.check_time_only``: Now also returns the signal. +- Bumped the minimum requirements to: + + - Python 3.11 + - NumPy 2.0 + - SciPy, Numba, libdlf, scooby (without minimum version) + v2.5.4: Bugfix ext. fcts with z+ up ----------------------------------- diff --git a/docs/manual/installation.rst b/docs/manual/installation.rst index bd97ac3..f95d958 100644 --- a/docs/manual/installation.rst +++ b/docs/manual/installation.rst @@ -13,8 +13,8 @@ or via ``pip``: pip install empymod -Requirements are the modules ``numpy``, ``scipy``, ``numba``, ``libdlf``, and -``scooby``. +Requirements are the modules ``numpy>=2.0.0``, ``scipy``, ``numba``, +``libdlf``, and ``scooby``. The modeller empymod comes with add-ons (``empymod.scripts``). These add-ons provide some very specific, additional functionalities. Some of these add-ons diff --git a/pyproject.toml b/pyproject.toml index cde7f73..df62918 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,13 +6,13 @@ build-backend = "setuptools.build_meta" name = "empymod" description = "Open-source full 3D electromagnetic modeller for 1D VTI media" readme = "README.rst" -requires-python = ">=3.10" +requires-python = ">=3.11" authors = [ {name = "The emsig community", email = "info@emsig.xyz"}, ] dependencies = [ - "numpy", - "scipy>=1.10", + "numpy>=2.0.0", + "scipy", "numba", "libdlf", "scooby", From 4444abdd23ce21d6290b42c739712c99df282781 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Sat, 31 Jan 2026 11:23:14 +0100 Subject: [PATCH 13/14] docs --- CHANGELOG.rst | 25 ++++++++++++----- docs/manual/info.rst | 67 ++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 79 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index b642c43..06fe95a 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -6,12 +6,16 @@ Version 2 ~~~~~~~~~ -v2.5.x +v2.6.x """""" -*latest* --------- +v2.6.0: Arbitrary waveforms +--------------------------- + +**2026-01-31** + +This version finally introduces the possibility for arbitrary waveforms. - Modelling routines: @@ -38,10 +42,13 @@ v2.5.x - New deprecations (will be removed in v3): - ``empymod.model.loop`` will be removed in v3. - - ``empymod.utils.check_bipole``: Now returns also the src/rec field and - type. - - ``empymod.utils.check_time``: Now also returns the signal. - - ``empymod.utils.check_time_only``: Now also returns the signal. + - ``empymod.utils.check_bipole``: New signature, it now returns also the + src/rec field and + type; old signature will be removed in v3. + - ``empymod.utils.check_time``: New signature, it now also returns the + signal; old signature will be removed in v3. + - ``empymod.utils.check_time_only``: New signature, it now also returns the + signal; old signature will be removed in v3. - Bumped the minimum requirements to: @@ -50,6 +57,10 @@ v2.5.x - SciPy, Numba, libdlf, scooby (without minimum version) +v2.5.x +"""""" + + v2.5.4: Bugfix ext. fcts with z+ up ----------------------------------- diff --git a/docs/manual/info.rst b/docs/manual/info.rst index bf2e917..35916ed 100644 --- a/docs/manual/info.rst +++ b/docs/manual/info.rst @@ -154,8 +154,8 @@ the example :ref:`sphx_glr_gallery_educational_dlf_standard_lagged_splined.py`. -Looping -------- +Looping to reduce memory +------------------------ By default, you can compute many offsets and many frequencies all in one go, vectorized (for the *DLF*), which is the default. The ``loop`` parameter gives @@ -218,6 +218,26 @@ A common alternative to this trick is to apply a lowpass filter to filter out the unstable high frequencies. +TEM / Loops +----------- + +The kernel of ``empymod`` computes the electric and magnetic Green's functions +for infinitesimal small dipoles. Combining various of these small dipoles +allows to model any kind of source and receiver configurations, including +loops, as frequently used in TEM (time-domain EM) measurements. To make the +correct assembly of dipoles is not always trivial, and the easiest way is +probably to look for a similar example and modify it. Here a list of loop +examples in the gallery: + +- :ref:`sphx_glr_gallery_fdomain_ipandq.py` (f-domain, but loop equipments, + GEM-2 & DUALEM-842); +- :ref:`sphx_glr_gallery_tdomain_ip_vrm.py`; +- :ref:`sphx_glr_gallery_tdomain_tem_temfast.py`; +- :ref:`sphx_glr_gallery_tdomain_tem_walktem.py`; +- :ref:`sphx_glr_gallery_educational_dipoles_and_loops.py`; +- :ref:`sphx_glr_gallery_reproducing_ward1988.py` (f- & t-demain). + + Hook for user-defined computation of :math:`\eta` and :math:`\zeta` ------------------------------------------------------------------- @@ -266,10 +286,45 @@ characteristics: And then you call ``empymod`` with ``res={'res': res-array, 'tau': tau, 'func_eta': my_new_eta}``. -Have a look at the corresponding example in the Gallery, where this hook is -exploited in the low-frequency range to use the Cole-Cole model for IP -computation. It could also be used in the high-frequency range to model -dielectricity. +Have a look at, e.g., the :ref:`sphx_glr_gallery_tdomain_cole_cole_ip.py` +example in the Gallery, where this hook is exploited in the low-frequency range +to use the Cole-Cole model for IP computation. It could also be used in the +high-frequency range to model dielectricity. + + +Hook for user-defined bandpass filter +------------------------------------- + +Similar to the hooks described above, there is since ``empymod v2.6.0`` a hook +to apply a bandpass filter (or any function) to the frequency domain result, by +proving a specific dictionary as ``bandpass`` argument. The signature of the +function must be ``func(inp, p_dict)``, where ``inp`` is the dictionary you +provide, and ``p_dict`` is a dictionary that contains all parameters so far +computed in empymod ``[locals()]``. Any change to the frequency domain result +must be done in-place, and the function does not return anything. Refer to the +time-domain loop examples in the gallery. The dictionary must contain at least +the keyword ``'func'``, containing the actual function, but can contain any +other parameters too. + +**Dummy example** + +.. code-block:: python + + def my_bandpass(inp, p_dict): + # Your computations, using the parameters you provided + # in `inp` and the parameters from empymod in `p_dict`. + # In the example line below, we provide, e.g., inp["cutoff"] + # Modification must happen in-place to "EM", nothing is returned. + + # Here just a simple cutoff to show the principle. + p_dict["EM"][p_inp["freq"] > inp["cutoff"] *= 0 + +And then you call ``empymod`` with ``bandpass={"cutoff": 1e6, "func": +my_bandpass}``. + +Have a look at, e.g., the :ref:`sphx_glr_gallery_tdomain_tem_walktem.py` +example in the Gallery, where this hook is applied. + Zero horizontal offset From e9dba2ebc56700f19a43551b8d7cb80d80413376 Mon Sep 17 00:00:00 2001 From: Dieter Werthmuller Date: Sat, 31 Jan 2026 11:28:07 +0100 Subject: [PATCH 14/14] docs II --- docs/manual/info.rst | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/docs/manual/info.rst b/docs/manual/info.rst index 35916ed..34c74c6 100644 --- a/docs/manual/info.rst +++ b/docs/manual/info.rst @@ -218,24 +218,22 @@ A common alternative to this trick is to apply a lowpass filter to filter out the unstable high frequencies. -TEM / Loops ------------ +Loops (FEM/TEM) +--------------- The kernel of ``empymod`` computes the electric and magnetic Green's functions for infinitesimal small dipoles. Combining various of these small dipoles allows to model any kind of source and receiver configurations, including -loops, as frequently used in TEM (time-domain EM) measurements. To make the -correct assembly of dipoles is not always trivial, and the easiest way is -probably to look for a similar example and modify it. Here a list of loop -examples in the gallery: +loops. To make the correct assembly of dipoles is not always trivial, and the +easiest way is probably to look for a similar example and modify it. Here a +list of loop examples in the gallery: -- :ref:`sphx_glr_gallery_fdomain_ipandq.py` (f-domain, but loop equipments, - GEM-2 & DUALEM-842); +- :ref:`sphx_glr_gallery_fdomain_ipandq.py` (GEM-2 & DUALEM-842); - :ref:`sphx_glr_gallery_tdomain_ip_vrm.py`; -- :ref:`sphx_glr_gallery_tdomain_tem_temfast.py`; -- :ref:`sphx_glr_gallery_tdomain_tem_walktem.py`; +- :ref:`sphx_glr_gallery_tdomain_tem_temfast.py` (TEM-FAST); +- :ref:`sphx_glr_gallery_tdomain_tem_walktem.py` (Walk-TEM); - :ref:`sphx_glr_gallery_educational_dipoles_and_loops.py`; -- :ref:`sphx_glr_gallery_reproducing_ward1988.py` (f- & t-demain). +- :ref:`sphx_glr_gallery_reproducing_ward1988.py`. Hook for user-defined computation of :math:`\eta` and :math:`\zeta`