bipole

empymod.model.bipole(src, rec, depth, res, freqtime, signal=None, aniso=None, epermH=None, epermV=None, mpermH=None, mpermV=None, msrc=False, srcpts=1, mrec=False, recpts=1, strength=0, **kwargs)

Return EM fields due to arbitrary rotated, finite length EM dipoles.

Calculate the electromagnetic frequency- or time-domain field due to arbitrary rotated, finite electric or magnetic dipole sources, measured by arbitrary rotated, finite electric or magnetic dipole receivers. By default, the electromagnetic response is normalized to source and receiver of 1 m length, and source strength of 1 A.

A magnetic dipole has a moment of \(I^m ds\). However, if the magnetic dipole is generated by an electric-wire loop, this changes to \(I^m = i\omega\mu A I^e\), where A is the area of the loop. The same factor \(i\omega\mu A\), applies to the receiver, if it consists of an electric-wire loop. This can be achieved by setting msrc or mrec to 'b'. The current implementation only handles loop sources and receivers in layers where \(\mu_r^h=\mu_r^v\); the horizontal magnetic permeability is used, and a warning is thrown if the vertical differs from the horizontal one.

Parameters:

src, reclist of floats or arrays

Source and receiver coordinates (m):

• [x0, x1, y0, y1, z0, z1] (dipole of finite length)

• [x, y, z, azimuth, dip] (dipole, infinitesimal small)

Note that when the coordinates are defined as finite length dipoles, it will still return an infinitesimal small dipole at its center unless srctpts;recpts is set to a value of 3 or higher.

For N sources or receivers, all variables must be of size N or 1 (in the latter case it will be expanded to N).

Angles (coordinate system is either left-handed with positive z down or right-handed with positive z up; East-North-Depth):

• azimuth (°): horizontal deviation from x-axis, anti-clockwise.

• +/-dip (°): vertical deviation from xy-plane down/up-wards.

Sources or receivers placed on a layer interface are considered in the upper layer.

depthlist

Absolute layer interfaces z (m); #depth = #res - 1 (excluding +/- infinity).

resarray_like

Horizontal resistivities rho_h (Ohm.m); #res = #depth + 1.

Alternatively, res can be a dictionary. See the main manual of empymod too see how to exploit this hook to re-calculate etaH, etaV, zetaH, and zetaV, which can be used to, for instance, use the Cole-Cole model for IP.

freqtimearray_like

Frequencies f (Hz) if signal==None, else times t (s); (f, t > 0).

signal{None, 0, 1, -1}, 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

anisoarray_like, default: ones

Anisotropies lambda = sqrt(rho_v/rho_h) (-); #aniso = #res.

epermH, epermVarray_like, default: ones

Relative horizontal/vertical electric permittivities epsilon_h/epsilon_v (-); #epermH = #epermV = #res. If epermH is provided but not epermV, isotropic behaviour is assumed.

mpermH, mpermVarray_like, default: ones

Relative horizontal/vertical magnetic permeabilities mu_h/mu_v (-); #mpermH = #mpermV = #res. If mpermH is provided but not mpermV, isotropic behaviour is assumed.

msrc, mrecbool or string, default: False

• True: If True, source/receiver (msrc/mrec) is magnetic, else electric.

• 'b': Source and receiver can also be set to msrc='b' and mrec='b', which means that they correspond to loops, simulating magnetic flux density B instead of magnetic field H, with B=μH. It computes the same as mrec=True, but multiplies it with μ of the layer of the receiver in the frequency domain (and accordingly four the source). Only implemented for isotropic magnetic permittivities at source/receiver level.

• 'j': The receiver can also be set to mrec='j', in which case the electric current density is returned (with the approximation that the current density is proportional to the electric field, J=σE). Only implemented for isotropic resistivities and electric permittivities at receiver level.

srcpts, recptsint, default: 1

Number of integration points for dipole source/receiver:

• srcpts/recpts < 3 : infinitesimal small dipole at the center.

• srcpts/recpts >= 3 : finite length approximation with N points.

strengthfloat, default: 0.0

Source strength (A):

• If 0, output is normalized to source and receiver of 1 m length (or 1 m² area if ‘b’), and source strength of 1 A.

• If != 0, output is returned for given source and receiver length, and source strength.

Note that the strength is simply a multiplication factor, and can be used for many other things therefore.

verb{0, 1, 2, 3, 4}, default: 2

Level of verbosity:

• 0: Print nothing.

• 1: Print warnings.

• 2: Print additional runtime and kernel calls

• 3: Print additional start/stop, condensed parameter information.

• 4: Print additional full parameter information

ht{‘dlf’, ‘qwe’, ‘quad’}, default: ‘dlf’

Flag to choose either the Digital Linear Filter (DLF) method, the Quadrature-With-Extrapolation (QWE), or a simple Quadrature (QUAD) for the Hankel transform.

htargdict, optional

Possible parameters depends on the value for ht:

• If ht='dlf':

  • dlf: string of filter name in empymod.filters or the filter method itself. (default: empymod.filters.Hankel().key_201_2009)

  • pts_per_dec: points per decade; (default: 0):

    • If 0: Standard DLF.

    • If < 0: Lagged Convolution DLF.

    • If > 0: Splined DLF

• If ht='qwe':

  • rtol: relative tolerance (default: 1e-12)

  • atol: absolute tolerance (default: 1e-30)

  • nquad: order of Gaussian quadrature (default: 51)

  • maxint: maximum number of partial integral intervals (default: 40)

  • pts_per_dec: points per decade; (default: 0)

    • If 0, no interpolation is used.

    • If > 0, interpolation is used.

  • diff_quad: criteria when to swap to QUAD (only relevant if pts_per_dec>0) (default: 100)

  • a: lower limit for QUAD (default: first interval from QWE)

  • b: upper limit for QUAD (default: last interval from QWE)

  • limit: limit for quad (default: maxint)

• If ht='quad':

  • rtol: relative tolerance (default: 1e-12)

  • atol: absolute tolerance (default: 1e-20)

  • limit: An upper bound on the number of subintervals used in the adaptive algorithm (default: 500)

  • a: Minimum wavenumber (default 1e-6)

  • b: Maximum wavenumber (default 0.1)

  • pts_per_dec: points per decade (default: 40)

ft{‘dlf’, ‘sin’, ‘cos’, ‘qwe’, ‘fftlog’, ‘fft’}, default: ‘dlf’

Only used if signal!=None. Flag to choose either the Digital Linear Filter method (Sine- or Cosine-Filter), the Quadrature-With-Extrapolation (QWE), the FFTLog, or the FFT for the Fourier transform. If ‘dlf’ it is ‘sin’ if signal>=0, else ‘cos’.

ftargdict, optional

Only used if signal!=None. Possible parameters depends on the value for ft:

• If ft='dlf', 'sin', or 'cos':

  • dlf: string of filter name in empymod.filters or the filter method itself. (Default: empymod.filters.Fourier().key_201_2012)

  • pts_per_dec: points per decade; (default: -1)

    • If 0: Standard DLF.

    • If < 0: Lagged Convolution DLF.

    • If > 0: Splined DLF

• If ft='qwe':

  • rtol: relative tolerance (default: 1e-8)

  • atol: absolute tolerance (default: 1e-20)

  • nquad: order of Gaussian quadrature (default: 21)

  • maxint: maximum number of partial integral intervals (default: 200)

  • pts_per_dec: points per decade (default: 20)

  • diff_quad: criteria when to swap to QUAD (default: 100)

  • a: lower limit for QUAD (default: first interval from QWE)

  • b: upper limit for QUAD (default: last interval from QWE)

  • limit: limit for quad (default: maxint)

• If ft='fftlog':

  • pts_per_dec: sampels per decade (default: 10)

  • add_dec: additional decades [left, right] (default: [-2, 1])

  • q: exponent of power law bias (default: 0); -1 <= q <= 1

• If ft='fft':

  • dfreq: Linear step-size of frequencies (default: 0.002)

  • nfreq: Number of frequencies (default: 2048)

  • ntot: Total number for FFT; difference between nfreq and ntot is padded with zeroes. This number is ideally a power of 2, e.g. 2048 or 4096 (default: nfreq).

  • pts_per_dec: points per decade (default: None)

  Padding can sometimes improve the result, not always. The default samples from 0.002 Hz - 4.096 Hz. If pts_per_dec is set to an integer, calculated frequencies are logarithmically spaced with the given number per decade, and then interpolated to yield the required frequencies for the FFT.

xdirectbool or None, default: False

Direct field calculation (only if src and rec are in the same layer):

• If True, direct field is calculated analytically in the frequency domain.

• If False, direct field is calculated in the wavenumber domain.

• If None, direct field is excluded from the calculation, and only reflected fields are returned (secondary field).

loop{None, ‘freq’, ‘off’}, default: None

Define if to calculate everything vectorized or if to loop over frequencies (‘freq’) or over offsets (‘off’). It always loops over frequencies if ht='qwe' or if pts_per_dec=-1. Calculating everything vectorized is fast for few offsets OR for few frequencies. However, if you calculate many frequencies for many offsets, it might be faster to loop over frequencies. Only comparing the different versions will yield the answer for your specific problem at hand!

squeezebool, default: True

If True, the output is squeezed. If False, the output will always be of ndim=3, (nfreqtime, nrec, nsrc).

Returns:

EMEMArray, (nfreqtime, nrec, nsrc)

Frequency- or time-domain EM field (depending on signal):

• If rec is electric, returns E [V/m].

• If rec is magnetic, returns H [A/m].

• If rec is b, returns B [T].

• If rec is j, returns J [A/m²].

EMArray is a subclassed ndarray with .pha and .amp attributes (only relevant for frequency-domain data).

The shape of EM is (nfreqtime, nrec, nsrc). However, single dimensions are removed, unless you set squeeze=False.

See also

dipole

EM fields due to infinitesimal small EM dipoles.

analytical

Analytical full- and half-space solutions.

Examples

In [1]: import empymod
   ...: import numpy as np
   ...: # x-directed dipole source: x0, x1, y0, y1, z0, z1
   ...: src = [-50, 50, 0, 0, 100, 100]
   ...: # x-directed dipole receiver-array: x, y, z, azimuth, dip
   ...: rec = [np.arange(1, 11)*500, np.zeros(10), 200, 0, 0]
   ...: # layer boundaries
   ...: depth = [0, 300, 1000, 1050]
   ...: # layer resistivities
   ...: res = [1e20, .3, 1, 50, 1]
   ...: # Frequency
   ...: freq = 1
   ...: # Calculate electric field due to an electric source at 1 Hz.
   ...: # [msrc = mrec = False (default)]
   ...: EMfield = empymod.bipole(src, rec, depth, res, freq, verb=3)
   ...: 

:: empymod START  ::  v2.5.5.dev12+g59de0f90f

   depth       [m] :  0 300 1000 1050
   res     [Ohm.m] :  1E+20 0.3 1 50 1
   aniso       [-] :  1 1 1 1 1
   epermH      [-] :  1 1 1 1 1
   epermV      [-] :  1 1 1 1 1
   mpermH      [-] :  1 1 1 1 1
   mpermV      [-] :  1 1 1 1 1
   direct field    :  Comp. in wavenumber domain
   frequency  [Hz] :  1
   Hankel          :  DLF (Fast Hankel Transform)
     > Filter      :  key_201_2009
     > DLF type    :  Standard
   Loop over       :  None (all vectorized)
   Source type     :  Electric field
   Receiver type   :  Electric field
   Source(s)       :  1 dipole(s)
     > intpts      :  1 (as dipole)
     > length  [m] :  100
     > strength[A] :  0
     > x_c     [m] :  0
     > y_c     [m] :  0
     > z_c     [m] :  100
     > azimuth [°] :  0
     > dip     [°] :  0
   Receiver(s)     :  10 dipole(s)
     > x       [m] :  500 - 5000 : 10  [min-max; #]
     > y       [m] :  0 - 0 : 10  [min-max; #]
     > z       [m] :  200
     > azimuth [°] :  0
     > dip     [°] :  0
   Required ab's   :  11

:: empymod END; runtime = 0:00:00.003079 :: 1 kernel call(s)


In [2]: EMfield[0]
Out[2]: np.complex128(1.6880934577857314e-10-3.083031298956568e-10j)