Document better loops

[prisae]

Any kind of loops (including central loops as a special case) are possible in empymod, but not trivial, and I get a lot of questions. Hence, it should be better documented.

Note: The loop function only converts a magnetic dipole measuring the magnetic field H [A/m] (theoretical, there are no magnetic dipole equipments) to a magnetic coil measuring the magnetic flux B [T].

There are examples of it:

• Figures 4.7 and 4.8 from Ward and Hohmann
• WalkTEM example
• In-phase and Quadrature, with GEM-2 and DUALEM-842 examples
• Difference between magnetic dipole and loop sources

For a central loop one needs to use bipole or dipole, loop does not work as source and receiver would be at the same location.

TODOs:

• Change in the above Ward and Hohmann example area/(radius/2) to circumference; it is confusing as it is now. Clarify Ward and Hohmann circular loop #233
• Create specifically a central loop example.
  • Simplest case: symmetric central loop
  • Square central loop
  • Central square loop within bigger square loop
  • Loop where the receiver is in an arbitrary place
  • Arbitrarily shaped loop
  • ...
  • Also document dipole H, coil B, and in time domain dB/dt
• Bring a TEMfast example from @aignerlukas, https://github.com/aignerlukas/TEMIP_sensitivity
• Maybe create a central loop function in empymod.model, as it is often requested.
• See comment below re point approximations of finite length dipoles. Maybe add to a gallery example.

[prisae]

import empymod
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

hoff = np.round(np.logspace(-1, 1, 6), 1) #[0.1, 0.5, 1, 2, 5]
inp = {
    'src': [0, 0, -0.5, +0.5, 0, 0],
    'rec': [hoff, np.zeros(len(hoff)), np.ones(len(hoff))*0.2, 0, 90],
    'depth': [],
    'res': 1.,
    'freqtime': 100,
    'mrec': True,
    'verb': 1,
    'xdirect': True,
}

points = [1, 3, 4, 5, 6, 7, 8, 9, 10, 1000]
res = np.zeros((len(points), len(hoff)))

for i, srcpts in enumerate(points):
    res[i, :] = empymod.bipole(**inp, srcpts=srcpts).amp()

error = 100*abs(res[:-1, :]-res[-1, :])/res[-1, :]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), constrained_layout=True)

ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.axhline(0, c='.6', zorder=1)
ax1.axvline(0, c='.6', zorder=1)
ax1.set_xlim([-1, max(hoff)+1])
ax1.set_ylim([-2, 2])
ax1.plot([0, 0], [-0.5, 0.5], c='k', lw=2, label='El. y-directed Source; z=0.0 m')
ax1.plot(hoff, np.zeros(len(hoff)), 'x', c='k', label='Mag. z-directed Receivers; z=0.2 m')
ax1.legend()

ax2.semilogy(error)
ax2.set_xticks(np.arange(len(points)-1), points[:-1])
ax2.set_xlabel('srcpts')
ax2.set_ylabel('Rel. Error (%)')

<img width="3600" height="1500" alt="Image" src="https://github.com/user-attachments/assets/4a6ef5be-23e7-4e44-b6aa-d5b89ec67589" />

ax2.legend(hoff, title='Offset (m)', loc=1)
ax2.axhline(1e0, c='.3', zorder=1)

fig.suptitle('Finite length dipole approximations with infinitesimall small (point) dipole sources\n'
             'Rule of thumb: if offset > 5x dipole length - point approximation is valid\n'
             'Rel. error calculated with a dipole of 1000 points (roughly mm segments).\n'
             'Analytical f-domain 1 Ωm fullspace solution for 100 Hz.')
fig.savefig('PointApproxOfFiniteLengthDipole', dpi=300)
empymod.Report()

[prisae]

LoopSketch.odp

[prisae]

An adaption of https://empymod.emsig.xyz/en/latest/gallery/educational/dipoles_and_loops.html

import empymod
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

def plot_result(data1, data2, x, title, vmin=-15., vmax=-7., rx=0):
    """Plot result."""
    #fig = plt.figure(figsize=(18, 10))
    fig = plt.figure(figsize=(15, 8))

    def setplot(name):
        """Plot settings"""
        plt.title(name)
        plt.xlim(rx.min(), rx.max())
        plt.ylim(rx.min(), rx.max())
        plt.axis("equal")

    # Plot Re(data)
    ax1 = plt.subplot(231)
    setplot(r"(a) |Re(srcpts=1)|")
    cf0 = plt.pcolormesh(rx, rx, np.log10(np.abs(data1.real)), linewidth=0,
                         rasterized=True, cmap="viridis", vmin=vmin, vmax=vmax,
                         shading='nearest')

    ax2 = plt.subplot(232)
    setplot(r"(b) |Re(srcpts=5)|")
    plt.pcolormesh(rx, rx, np.log10(np.abs(data2.real)), linewidth=0,
                   rasterized=True, cmap="viridis", vmin=vmin, vmax=vmax,
                   shading='nearest')

    ax3 = plt.subplot(233)
    setplot(r"(c) Error real part")
    error_r = abs(data1.real-data2.real)/(abs(data1.real)+abs(data2.real))*200
    cf2 = plt.pcolormesh(rx, rx, np.log10(error_r), vmin=-2, vmax=2,
                         linewidth=0, rasterized=True,
                         cmap=plt.get_cmap("RdBu_r", 8),
                         shading='nearest')

    # Plot Im(data)
    ax4 = plt.subplot(234)
    setplot(r"(d) |Im(srcpts=1)|")
    plt.pcolormesh(rx, rx, np.log10(np.abs(data1.imag)), linewidth=0,
                   rasterized=True, cmap="viridis", vmin=vmin, vmax=vmax,
                   shading='nearest')

    ax5 = plt.subplot(235)
    setplot(r"(e) |Im(srcpts=5)|")
    plt.pcolormesh(rx, rx, np.log10(np.abs(data2.imag)), linewidth=0,
                   rasterized=True, cmap="viridis", vmin=vmin, vmax=vmax,
                   shading='nearest')

    ax6 = plt.subplot(236)
    setplot(r"(f) Error imag part")
    error_i = abs(data1.imag-data2.imag)/(abs(data1.imag)+abs(data2.imag))*200
    plt.pcolormesh(rx, rx, np.log10(error_i), vmin=-2, vmax=2,
                   linewidth=0, rasterized=True,
                   cmap=plt.get_cmap("RdBu_r", 8),
                   shading='nearest')

    # Colorbars
    fig.colorbar(cf0, ax=[ax1, ax2, ax3], label=r"$\log_{10}$ Amplitude (A/m)")
    cbar = fig.colorbar(cf2, ax=[ax4, ax5, ax6], label=r"Relative Error")
    cbar.set_ticks([-2, -1, 0, 1, 2])
    cbar.ax.set_yticklabels([r"$0.01\,\%$", r"$0.1\,\%$", r"$1\,\%$",
                             r"$10\,\%$", r"$100\,\%$"])

    # Axis label
    fig.text(0.4, 0.05, "Inline Offset (m)", fontsize=14)
    fig.text(0.08, 0.5, 'Crossline Offset (m)', rotation=90, fontsize=14)

    # Title
    fig.suptitle(title, y=.95, fontsize=20)
    return fig

zsrc = 1.31

# Model
model = {
    'depth' : [0., 5., 7.],
    'res' : [1e-14, 1., 0.2, 1.],
    'freqtime': 1000,    # 100 Hz
    'htarg': {'pts_per_dec': -1},
    'verb': 1,
}

# Survey parameters

for factor in [10, 100]:

    x = ((np.arange(502))-250.5)/factor
    rx = np.repeat([x, ], np.size(x), axis=0)
    ry = rx.transpose()
    rxx = rx.ravel()
    ryy = ry.ravel()
    
    rec_dip = [rxx, ryy, 1.21]
    
    square_pts = +empymod.dipole(src=[+0.5, +0.0, zsrc], rec=rec_dip, ab=62, **model).reshape(np.shape(rx))
    square_pts -= empymod.dipole(src=[+0.0, +0.5, zsrc], rec=rec_dip, ab=61, **model).reshape(np.shape(rx))
    square_pts -= empymod.dipole(src=[-0.5, +0.0, zsrc], rec=rec_dip, ab=62, **model).reshape(np.shape(rx))
    square_pts += empymod.dipole(src=[+0.0, -0.5, zsrc], rec=rec_dip, ab=61, **model).reshape(np.shape(rx))
    
    inp_dip = {
        'rec': [rxx, ryy, 1.21, 0., 90.],
        'mrec': True,
        'srcpts': 5  # Gaussian quadr. with 5 pts to simulate a finite length dip.
    }
    
    square_dip = +empymod.bipole(src=[+0.5, +0.5, -0.5, +0.5, zsrc, zsrc], **inp_dip, **model)
    square_dip += empymod.bipole(src=[+0.5, -0.5, +0.5, +0.5, zsrc, zsrc], **inp_dip, **model)
    square_dip += empymod.bipole(src=[-0.5, -0.5, +0.5, -0.5, zsrc, zsrc], **inp_dip, **model)
    square_dip += empymod.bipole(src=[-0.5, +0.5, -0.5, -0.5, zsrc, zsrc], **inp_dip, **model)
    square_dip = square_dip.reshape(np.shape(rx))
    
    fig = plot_result(square_pts, square_dip, x, 'Four points vs four finite dipoles', vmin=-8, vmax=1, rx=x)